text:
# The Liverpool-Edinburgh High Proper Motion Survey

## 1 Introduction

### 1.1 Proper Motion Surveys

Proper motion surveys can give us a good insight into stellar populations and dynamics in the Solar neighbourhood. They identify nearby objects, allow analysis of their space motions and, with data in more than one passband, give population differentiation. These three factors mean that proper motion surveys are ideal in addressing the problems of the local stellar mass function and the halo mass function, both of which are currently poorly constrained (e.g. Henry et al., 1997, Lee 1993). The problem with the local stellar mass function is in large part due to the undersampling of nearby stars upon which it is mainly based. Identification of objects such as the nearby brown dwarf Kelu-1 (Ruiz et al. 1997), which was discovered as part of the Calan-ESO proper motion survey (Ruiz et al. 2001) and has an astrometrically estimated distance of $``$10pc, are vital in this respect. Both present day and initial halo mass functions are poorly constrained (Lee 1993) due to a number of factors. These include the relatively low number of nearby halo subdwarf stars and the difficulty in obtaining accurate distances for them. However, halo stars have greater space velocity and velocity dispersion than disc stars which means that a halo star will have a higher proper motion than a disc star at the same distance. When combined with their subluminosity this provides an excellent method of differentiating them from disc stars (seeSect. 5$``$1).

However, all of the existing surveys suffer from one or more of the following problems;

* internal inhomogeneities which make statistical analysis difficult and unreliable,
* low areal coverage which makes a survey subject to low number statistics,
* use of non-contiguous survey fields which prevents cross-checking of discovered objects and completeness levels from areas where contiguous fields overlap.

One of the main purposes of this survey is to address these problems by conducting an internally homogeneous large area survey in contiguous fields by means of scripted computer algorithms. By using only the R-band plates we also ensured that the survey would have a high completeness relative to other surveys which require an object to be present in more than one passband, thus missing faint objects with extreme colours, such as brown dwarfs and ancient halo subdwarfs and white dwarfs. This is important because these objects are among the most astrophysically interesting to be selected by proper motion surveys. Using only the R-band plates also means that our completeness is easier to assess.

### 1.2 The Luyten Catalogues

The largest catalogue of high proper motion stars to date is the New Luyten Catalogue of Stars with Proper Motions Greater Than Two Tenths of an Arcsecond (NLTT) (Luyten, 1979). The NLTT and the Luyten Half Second Proper Motion Catalogue (LHS) (Luyten, 1979) were compiled over a period of $``$60 years from two main surveys; the Bruce Proper Motion Survey (BPM), covering an area from the South Celestial Pole to a declination of $``$+25, and the Palomar Observatory Sky Surveys (POSS I and II) covering the area from the North Celestial Pole to $``$–33. The nominal faint limit of the Bruce survey is R$``$16$``$5mag and the faint limit of the POSS is R$``$20mag (Morgan et al., 1992), the Luyten Catalogues have no stated bright limit. The Luyten proper motion measurements were made using two different techniques; (i) manual plate blinking and (ii) automatic machine detection (Williams, 2000). The Luyten catalogues therefore contain some internal inhomogeneities such as a brighter limiting magnitude for fields south of –33 and variations in the completeness of the catalogues between different areas of the surveys. The stated proper motion limits of the NLTT are 0$``$18$`\mathrm{}`$yr-1$`\mu `$2$``$5$`\mathrm{}`$yr-1 whilst the LHS comprises stars from the NLTT with $`\mu `$0$``$5$`\mathrm{}`$yr-1.

### 1.3 Other Surveys and the Completeness of the Luyten Catalogues

There have been a number of estimates of the completeness of the NLTT and LHS catalogues in recent years, mostly based on the results of small area proper motion surveys.

Scholz et al. (2000) conducted a survey in 40 fields using the Automatic Plate Measuring microdensitometer (APM) at Cambridge, UK. They discuss the completeness of previous catalogues, as well as how they find far more stars in the relatively incomplete southern area of the NLTT, but beyond stating that the NLTT is incomplete they do not quantify their own completeness or, using their data, the completeness of the NLTT.

Evans (1992) compared the results from four fields of the APM proper motion project to the results for the same areas from the NLTT. He discovered a number of NLTT stars which, when measured with the APM microdensitometer, proved to have proper motions below the 0$``$18$`\mathrm{}`$yr-1 cutoff limit. He also found a larger number of stars which had proper motions above this limit which were not listed in the NLTT. By conducting a Monte Carlo simulation he showed that these discrepancies could be almost entirely attributed to random measuring errors. He further stated that these errors led to a contamination of between 10$`\%`$ and 20$`\%`$ and an incompleteness of $``$16$`\%`$ (ie a completeness of $``$84$`\%`$). However, whilst the comparison was meticulously carried out and extremely detailed it relied on comparison of a very small number of fields. It is therefore possible that it suffers from selection effects and small number statistics.

Ruiz et al. (2001) compared their results from the Calan-ESO proper motion catalogue to the Luyten catalogues and concluded that for $`\mu `$$``$0.8$`\mathrm{}`$yr-1 and R$``$13mag the LHS catalogue is $``$60$`\%`$ complete. They also state that the Luyten Two Tenths catalogue (LTT) is $``$40$`\%`$ complete for R$``$13mag, but this refers to the 1975 version of the Two Tenths catalogue, which was superseded by the more complete NLTT in 1979. The Ruiz et al. comparison suffers from a similar problem to that of Evans in that it used the relatively small number of 14 fields. In addition to this their survey was conducted by manual blink comparison which is always prone to human error, as highlighted by the larger incompleteness of the southern BPM section of the NLTT over the northern POSS section. It was also conducted at low Galactic latitude $``$b$``$$``$40 because fields at high Galactic latitude do not contain enough background stars for efficient blink comparison, which could be an important contributory factor in the incompleteness of the BPM section of the NLTT.

Dawson (1986) applied the V/Vmax method as described by Schmidt (1968, 1975) to estimate a completeness of $``$90$`\%`$ for $`\delta `$$``$-33, $``$b$``$$`>`$10, R$``$18mag and $`\mu `$$``$0$``$5$`\mathrm{}`$yr-1 (i.e. the northern areas of the LHS).

Flynn et al. (2001) applied a statistical test to the NLTT consisting of taking two concentric spheres centred around the Sun whose volumes are in the ratio 2:1 (thereby giving them a radial ratio of 1$``$259 which corresponds to an average magnitude difference of 0$``$5mag between the shells) and assuming that the completeness for stars in the inner sphere is 100$`\%`$. It is then possible to calculate the completeness of the stars in the outer ‘shell’, and by applying this to different magnitude bins they derive a plot of completeness for the LHS and NLTT catalogues relative to the completeness at R=13 (assumed for the purposes of the plot to be 100$`\%`$ complete). This analysis shows the completeness of the NLTT dropping steadily from 80$`\%`$ at R=14 to 60$`\%`$ at R=18$``$5mag. However, Monet et al. (2000) suggest that this analysis is flawed since it does not take into account the different space densities of stars at high and low Galactic Latitude. Their recomputation of the Flynn et al. test gives a completeness at low Galactic latitude (15$``$$``$b$``$$``$35) of between 75$`\%`$ and 85$`\%`$ for 14$``$R$``$18mag.

Monet et al. (2000) conducted a survey in 35 fields using POSS II plates and estimated that the LHS is $``$90$`\%`$ complete, as well as recalculating the result of Flynn et al. (discussed above).

## 2 The Data and Search Process

### 2.1 The SuperCOSMOS Sky Survey

The SuperCOSMOS Sky Survey (SSS) (Hambly et al. 2001a) is a project of the Wide-Field Astronomy Unit at the Royal Observatory Edinburgh. The purpose of the project is to digitise Schmidt survey plates from the ESO, UK and Palomar (southern area) Schmidt surveys with the SuperCOSMOS microdensitometer machine and make the digitised data available online to the general astronomical community.

### 2.2 The Survey Data

The proper motion survey was conducted by comparison of two epochs of SuperCOSMOS digitized R-band Schmidt plates. The plates used are the ESO Schmidt and the Anglo-Australian UK Schmidt R-band surveys. The two surveys were taken at different observatories by different groups which means that the epoch differences between the plates vary. In the majority of cases the ESO plate was taken before the UK plate, although in a very small number of cases the UK plate is earlier. The average epoch difference is 8$``$5 years and the vast majority of fields have epoch differences between 5 and 13 years (see Fig. 2). The survey was conducted using only the two R-band plates in order not to miss faint objects with extreme colours, such as late M dwarfs or brown dwarfs, which might not be found by a search which also required detections in the BJ and I band plates. Once the survey was completed the positions of the candidate high proper motion stars were extrapolated to the BJ and I band plates in order to extract BRI colour information where possible.

### 2.3 The Digitized Data

The data has been digitized using the SuperCOSMOS microdensitometer at the Royal Observatory Edinburgh UK (Hambly et al. 1998) as part of the SuperCOSMOS Sky Survey (Hambly et al., 2001, a, b & c). Each object on the plate is reduced to 32 image parameters which describe its position, size, intensity, etc., in a format known as Image Analysis Mode (IAM) (Hambly et al., 2001a). These parameters make it possible to perform full astrometry and photometry on the data and to compare the data from two or more plates of the same field.

### 2.4 The Search Programs

The survey was conducted by a series of programs which manipulate the IAM files to perform full relative astrometry on the two epochs of each field. The programs attempt to pair all objects between the two plates from which the high proper motion stars may be selected. Objects are selected for the search process on the basis of a number of the IAM image parameters;

* the quality flag (which indicates if there are any problems with the object such as it being too faint, too bright, too large or too small),
* the deblending flag which indicates if the object is part of a larger object (e.g. part of the diffraction spike or halo of a bright star),
* the classification flag, which indicates whether the object is a galaxy, a star, dust or dirt, or too faint to classify.

The first task is to match the coordinate system of the UK plate to that of the ESO plate and trim the edges from the data (the latter being important because the UK plates are 6 square whilst the ESO plates are 5 square). The data for both plates is then reordered within the IAM files into 24 lanes so that the pairing programs only need to search a limited area of the plates rather than the whole of the dataset thus saving computer time and resources.

The position of each object on the UK plate is then checked against the ESO plate for any object at the same position. This is a simple first pass procedure to identify objects which have $``$ zero proper motion. Every ESO plate object left unpaired by this procedure is searched for on the UK plate in annuli of gradually increasing radius. The search looks for any object which has similar IAM image parameters to that of the object on the ESO plate. When a match is found an identifier is written to a file of paired objects and the program moves on to the next object. If no match is found within the cutoff radius the program writes a null identifier to the pairs file.

The search is conducted in this two step manner because the likelihood of finding a match decreases as the search radius increases. Searching the entire area out to the cut-off radius in one pass would use a large amount of computer time and processing power, which would be mostly wasted since for the vast majority of stars a match will be found within the first few annuli. Additionally, searching the full radius in one go could result in two or three objects being identified as possible pairings for every ESO object, requiring the software to choose between these objects as to which is the most likely to be the correct one. Again, this would again be a waste of time and processor power, since the most likely match of two or more similar candidates would always be the one nearest to the first epoch position.

We then perform full relative astrometry on the stars in x and y plate coordinates (measured in mm). The next stage is to calculate a transformation from x-y to RA-dec coordinates, which can then be applied to the x-y relative astrometry to yield full relative proper motions in RA and dec. Once this has been done catalogue stars can be selected by applying proper motion cuts to the paired objects. Finally the catalogue can be compared to an electronic version of the NLTT (see Sect. 4).

Limits of 0$``$18$`\mathrm{}`$yr-1$`\mu `$20$``$0$`\mathrm{}`$yr-1, and 9$``$R$``$19$``$5mag were applied to the final selection of proper motion objects. The lower proper motion limit was chosen to avoid any erroneous pairings due to the positional errors in the SuperCOSMOS data and to be consistent with the NLTT. The high proper motion limit was chosen to be well above any expected proper motion in order to include all potential objects. The bright magnitude limit was chosen because of the problems encountered with bright objects on a Schmidt plate (e.g. size of image, diffraction and saturation), whilst the faint magnitude limit was chosen because it becomes extremely difficult to differentiate by eye between dust and true stellar images for R$`>`$19$``$5mag on a Schmidt plate, which must always be the final test for any proper motion object discovered from the digitized data. This has the additional effect that the object classification contained in the IAM data becomes less accurate for faint objects which also makes selection of good pairs by use of the programs extremely difficult for images fainter than R=19$``$5mag.

After the proper motion selection was completed the data was compared to the BJ and I band plates (where available) in order to obtain magnitude data for the objects in those passbands.

## 3 The Catalogue

The catalogue comprises 6605 objects, 399 of which are extra detections of stars from the plate overlap regions and are already listed in the catalogue. This means that there are actually 6206 high proper motion stars in the catalogue. They lie within an area of -82.5$``$dec$``$-22.5 and 21hrs$``$RA$``$5hrs with magnitude limits of 9$``$mR$``$19.5mag and proper motion limits of 0.18$``$$`\mu `$$``$20$`\mathrm{}`$yr-1. Table 1 is a sample from the catalogue, comprised of the first 40 objects.

The columns are; catalogue number, right ascension (hr, m, s) and declination (, $`\mathrm{}`$, $`\mathrm{}`$) (equinox J2000 at the epoch of the ESO plate), multiple detection flag, ESO/SRC number of the field the object was located on, Epoch of ESO observation, proper motion in R.A., error on proper motion in R.A., proper motion in dec, error on proper motion in dec, total proper motion, error on total proper motion (all proper motions and their associated errors are given in $`\mathrm{}`$yr-1), position angle of proper motion (in degrees), apparent magnitudes: BJ, RESO, RUK, I. The multiple detection flag indicates where, in an overlap region between plates, an object is identified more than once (in some cases four times). This flag is 0 if the object is not a multiple detection or has the more accurate data of the mutiple detections. The flag is 1 if the object is a multiple detection with less accurate data, e.g. where there is no matching object was identified on the BJ or I band plates or the proper motion error is higher than for other detections. In those cases where no match could be found in the BJ or I band plates the apparent magnitude is given as 0.00.

## 4 Completeness of the Survey

### 4.1 The completeness test

Our completeness test consists of placing 100 simulated proper motion stars, each of which has the same magnitude and proper motion, into the data set for an ESO field and conducting the search process to see how many are recovered. This is repeated for the same magnitude with increasing proper motion starting from 0$``$3$`\mathrm{}`$yr-1 then 0$``$5$`\mathrm{}`$yr-1 and then in steps of 0$``$5$`\mathrm{}`$yr-1 up to 10$`\mathrm{}`$yr-1. The same procedure is then repeated for every 0$``$5 magnitudes from R=12mag to R=19mag.

The lower limit of 0$``$3$`\mathrm{}`$yr-1 was used because the simulated stars are added to the data after the coordinate matching procedure but before the astrometric solution is applied to the data. The calculated proper motions of the simulated stars can therefore vary by more than 10$`\%`$, which means that an input proper motion of 0$``$2$`\mathrm{}`$yr-1 results in some of the simulated stars being lost due to their calculated proper motions falling below the lower proper motion cutoff limit of 0$``$18$`\mathrm{}`$yr-1. The upper limit of 10$`\mathrm{}`$yr-1 was chosen because no candidate star was discovered by the survey with a proper motion above $``$10$`\mathrm{}`$yr-1.

One of the main problems with assessing the overall completeness of this survey is the way in which the plate epoch differences vary from field to field. An increase in epoch difference means that, for a given proper motion, the physical image separation of a star between the two epochs will also increase. This results in an increase in the probability of an erroneous pairing which will in turn reduce the completeness. Fortunately this problem can be quantified by running the test on fields with different epoch differences. For this reason the tests were carried out on three fields;

* Field 188, with a 2$``$980 year epoch difference,
* Field 479, with a 9$``$078 year epoch difference,
* Field 158, with a 16$``$911 year epoch difference,

These fields are representative of the shortest and longest epoch differences as well as an epoch difference roughly halfway between these two extremes.

### 4.2 Epoch difference limits on fields used in the survey

Fields with epoch differences of less than $``$3 years were discarded due to a large increase in false pairings at these shorter epoch differences. It was originally thought that these false pairings were caused by dirt and dust which are randomly distributed on the plates. However, there are too few dust particles on the UK plates which could be paired with the stars on the ESO plates to have such a large effect on the number of pairings. The overall proper motion errors are, on average, about 0.04$`\mathrm{}`$yr-1, but analysis of these errors reveals that they increase for fields with a short epoch difference. The proper motion error is propogated from the centroiding errors on the SuperCOSMOS data. This increase begins to have a significant effect on the proper motions at an epoch difference of slightly less than 3 years. The falsely identified stars are therefore not wrongly paired but are the result of the higher proper motion errors on stars with proper motions which are actually below the 0.18$`\mathrm{}`$yr-1 cut off limit in plates with short epoch differences.

Although the completeness decreases for plates with large epoch differences (see Fig. 3) this is simply a systematic effect of the increase in the plate to plate separations of stars at any given proper motion increasing the likelihood of false pairings. Since this only affects faint objects with relatively large proper motions no upper limit was applied to the allowed epoch difference of plate pairs.

### 4.3 Completeness test result

Fig. 3 shows that for fields with a short epoch difference the test completeness remains over 90% for all magnitudes and proper motions tested. Fig. 3 shows how the completeness decreases for increasing epoch difference and Fig. 4 how it decreases for R$`>`$18.5mag. This drop off is due partly to the increasing number of stars at faint magnitudes, and partly to dust and dirt on the plates, which can be difficult to distinguish from faint stars. The former results in an increase in the density of faint stars which in turn increases the probability of a false pairing. The latter is only a small effect for the reasons discussed above. Morgan et al. (1992) give the limiting magnitudes of both ESO R and UK R plates as 22, although they point out that the limiting magnitudes of Schmidt plates are not defined in a consistent manner. Hambly et al. (2001b) give the limit for identifying an object above the noise of an ESO Schmidt plate as R$``$21. Our limit of R=19.5mag is therefore a little conservative, but appropriate since we need to do more than just identify objects above the noise and are aiming for a uniform catalogue.

### 4.4 Cumulative number of proper motion stars

Another measure of the completeness of the survey is relationship between the proper motion and the cumulative number of stars with proper motions equal to or greater than a given proper motion. For any tangential velocity, VT,

$$\mu r^1$$
(1)
the volume (V) is proportional to r3.
For a uniform distribution of stars the number of stars per unit volume is constant (i.e. $`\rho `$ is constant). Since

$$Nc=\rho V$$
(2)
this implies that

$$Ncr^3$$
(3)
and, therefore, that

$$Nc\mu ^3.$$
(4)
Fig. 5 is plot of Nc against $`\mu `$ for all stars in this survey with $`\mu `$$``$0$``$2$`\mathrm{}`$yr-1. The data has a slope of -3.007 for 0$``$2$``$$`\mu `$$``$1$``$0$`\mathrm{}`$yr-1 and demonstrates that this survey is complete within these limits.

### 4.5 Completeness of the NLTT

This survey identified 6206 high proper motion objects in $``$7$``$5$`\%`$ of the whole sky with a completeness of $``$90$`\%`$. This implies $``$86,400 stars with magnitudes 9$``$R$``$19.5mag and proper motions $`\mu `$$``$0.18$`\mathrm{}`$yr-1 over the whole sky compared to the 58,845 stars in the NLTT, which would suggest at first sight that the overall completeness of the NLTT is $``$70$`\%`$. This is, however, extremely misleading since the NLTT is made up of two surveys (see Sect. 1$``$2). The POSS plate section of the survey (north of dec$``$-33) has a nominal faint magnitude limit of R$``$18$``$5mag, and was mostly conducted by a laser microdensitometer machine, thus making it largely free from human error. The BPM section (south of $``$-33) has a nominal faint R-band magnitude limit of R$``$16$``$5mag and was conducted entirely by eye on a blink comparator (which Luyten himself developed). The Luyten catalogues are therefore far less complete in the southern BPM area, as well as having a brighter limiting magnitude. Figs. 6 and 7 show how the completeness and faint limits of the different sections of the NLTT vary. Fig. 6 shows the number of high proper motion stars in the NLTT and this survey for the 43 fields in this survey which have field centres with dec$``$\- 30. Fig. 7 shows the number of high proper motion stars in the NLTT and this survey for the 73 fields in this survey which have field centres with dec$``$-40. The overlap between fields with centres at dec=-35 and dec=-30 means that the fields with centres at dec=-35 contain a number of stars discovered from the POSS plates which would increase the number of stars in the NLTT relative to this catalogue. This would give a false indication of the relative numbers of stars between the two catalogues, and for this reason the 15 fields with plate centres at dec=-35 were omitted from this plot.

As discussed in Sect. 1$``$3, Monet et al. have also shown that the completeness of the NLTT is dependent on Galactic latitude. Our survey was conducted at Galactic latitudes -30$``$b$``$-90 which most nearly corresponds to the high Galactic latitude fields which Monet et al. show to be less complete. Our analysis shows that in the POSS section of the NLTT is $``$85$`\%`$ complete for 10$``$5$``$R$``$19$``$5mag and that the BPM section is $``$90$`\%`$ complete for R$``$13mag with the completeness dropping off rapidly for fainter magnitudes.

## 5 Object populations and selection for follow up observations

As noted in Sect. 1$``$1, proper motion surveys can give an insight into both the local stellar neighbourhood and halo populations. In order to have a high proper motion an object must be either nearby or have a high space velocity (or both). High proper motion surveys therefore naturally select nearby and halo stars. Analysis of the objects identified by the survey requires differentiation of the populations, and this can be achieved by use of reduced proper motion and colour-colour diagrams. Astrophysically interesting objects can then be selected for follow up observations on this basis.

We recently obtained data for a number of stars selected from the catalogue on the ESO New Technologies Telescope at the La Silla observatory, Chile. The stars were selected as; those with the highest proper motions, candidate brown dwarfs, halo subdwarfs. These data are currently being reduced.

### 5.1 Reduced Proper Motion

The concept of reduced proper motion (H) was first used by Ejnar Hertzsprung (Luyten, 1968) as an alternative to absolute magnitude in order to produce ‘statistical’ HR diagrams. It combines the observable properties of proper motion, $`\mu `$, and apparent magnitude, $`m`$, via the equation,

$$H=m+5+5log\mu $$
(5)
From the definition of parallax it can be seen that,

$$V_T=4.74\frac{\mu }{\pi }$$
(6)
where $`V_T`$ is the tangential velocity measured in kms-1, and $`\pi `$ is the parallax in arcseconds, and since,

$$\pi =\frac{1}{d}$$
(7)
where $`d`$ is in parsecs, it follows that,

$$d=\frac{V_T}{4.74\mu }.$$
(8)
Using this and the relation between apparent and absolute magnitudes,

$$m=M+5log\left(\frac{d}{10}\right),$$
(9)
we can substitute into eqn.5 to get,

$$H=M+5logV_T3.379.$$
(10)
From this it can be seen that the reduced proper motion, as calculated from the measured properties of proper motion and apparent magnitude, is equivalent to a combination of the absolute magnitude and tangential velocity, which are intrinsic properties of a star. Plotting the reduced proper motion against colour thus gives a plot which is in essence an H-R diagram modified to include tangential velocity, called a reduced proper motion diagram.

One of the advantages of such a plot is that the higher mean tangential velocity of halo subdwarfs, when combined with their intrinsic subluminosity makes them far easier to differentiate from the disk population. However, the spread in H of objects in a reduced proper motion diagram is far larger than the spread in magnitude in a standard H-R diagram because of the random spread of velocities about the mean values for the different populations. It also means that stars can be misidentified due to unusually high or low tangential velocities.

Fig. 8 is a plot of the B-band reduced proper motion (HB) against B-I and shows the differentiation between the white dwarf, halo and disc populations. It should be noted that for $``$150 faint red stars no matches could be found on the BJ plates.

### 5.2 The Colour-Colour Diagram

The BRI two colour diagram (Fig. 9) utilises the BJ and I plates from the UK Schmidt survey. Halo subdwarfs and giant stars occupy the bottom right of this plot, the white dwarfs form the tail to the lower left and the disc main sequence runs from above the halo subdwarf region to the top right of the plot. This plot can therefore be used in conjunction with the reduced proper motion plots to differentiate between the various populations of objects.

## 6 Discussion

Over the past few years there have been several papers which offer a variety of values for the completeness of the Luyten catalogues, and discuss its shortcomings (see Sect. 1$``$2). The values offered for the completeness of the Luyten catalogues vary from 60$`\%`$ for mR$`>`$13 (Ruiz et al., 2001,comparison to Calan-ESO survey conducted at low galactic latitude) to 90$`\%`$ for mR$``$18 and $`\mu `$$``$0$``$5$`\mathrm{}`$yr-1 (Dawson 1986, V/Vmax method applied to northern areas of LHS)

Our completeness levels are the first to be calculated from an analysis of the search technique itself, and show that for the vast majority of stars in the catalogue (R$``$19mag and 0$``$2$``$$`\mu `$$``$2$``$5$`\mathrm{}`$yr-1) the completeness is greater than 90$`\%`$. Fig. 6 therefore indicates that for the Northern areas ($`\delta `$$`>`$-33, $``$b$``$$`>`$10) the NLTT is $``$85$`\%`$ complete. Fig. 7 shows the incompleteness of the Southern part of the NLTT for R$`>`$13mag ($``$25$`\%`$ complete for 15$``$R$``$16$``$5mag and $``$15$`\%`$ for 16$``$5$``$R$``$18mag).

We used the HB, B-I and the B–R, R–I plots to select halo subdwarfs. This selection yielded 1189 halo subdwarfs from a total of 6104 stars which have both B and I band magnitudes ($``$19.5$`\%`$). Table 6 of Reid (1984) gives four model predictions of the relative numbers of disc, halo, intermediate (thick disc) and white dwarf populations which would be discovered by a polar proper-motion survey. The models are based on different combinations of the disc and halo luminosity functions and the kinematics of the halo. Our result gives a value which falls between the predictions of models A and D and excludes models B amd C. Model A is based largely on the results of Wielen (1974), whilst model D uses the Luyten (1968) disc luminosity function, a ‘mean globular’ halo luminosity function and the halo kinematics of Oort (1965). This analysis is described in more detail in Pokorny, Jones & Hambly, 2002 (in preparation).

## Acknowledgments

The authors would like to thank all of the staff of the SuperCOSMOS project at the ROE and Mike Read for producing Fig. 1. RSP would like to thank J. M. Porter, P. A. James and A. M. Newsam for helpful discussions and suggestions. RSP acknowledeges support of an LJM University fellowship.


format: text

text:
# Raft Instability of Biopolymer Gels*footnote **footnote * Published in Phys. Rev. Lett. 87, 158101 (2001)

## Abstract

Following recent X-ray diffraction experiments by Wong, Li, and Safinya on biopolymer gels, we apply Onsager excluded volume theory to a nematic mixture of rigid rods and strong “$`\pi /2`$” cross-linkers obtaining a long-ranged, highly anisotropic depletion attraction between the linkers. This attraction leads to breakdown of the percolation theory for this class of gels, to breakdown of Onsager’s second-order virial method, and to formation of heterogeneities in the form of raft-like ribbons.

The demonstration by Lars Onsager in 1949  that solutions of long thin rods undergo a first-order phase-transition from an isotropic to a birefringent nematic phase has remained a landmark achievement of theoretical statistical physics. He showed that at the isotropic/nematic transition point the rod volume fraction $`\varphi `$ is surprisingly low, of the order of the aspect ratio $`D/L`$ of the rods, $`D`$ being the rod diameter and $`L`$ the rod length. Stiff biopolymers provide interesting applications of Onsager theory. In particular, the rod-like filamentous protein actin, which carries many biophysical functions , exhibits an isotropic-to-nematic transition  at a critical volume fraction in approximate agreement with the Onsager criterion . Inside cells, micron-size actin filaments are part of a gel structure, the cytoskeleton, with the filaments cross-linked by reversibly bound proteins like $`\alpha `$-actinin. The unusual elastic properties of in-vitro actin gels  have recently been the focus of an intense theoretical effort .

In-vitro studies of the sol-gel transition of fixed-length actin filaments are currently interpreted in terms of percolation theory , in which linkers are assigned, at random, to rods providing connections to neighboring rods with no preferred crossing angle. The sol-gel point, where the real and imaginary parts of the elastic moduli are comparable in magnitude, is identified as the percolation threshold, i.e., the linker concentration for which a connected “path” of linked rods stretches across the sample for the first time. An important validity condition is that the concentration of “native” rod-rod contacts - i.e., the mean number density of rod-rod contacts before linkers are added - exceeds the percolation threshold, otherwise the reduction in configurational freedom of the rods imposed by the linkers would lead to phase-separation. The number of native contacts per rod in a dilute isotropic solution of rods is of the order of $`(L/D)\varphi `$ so a rod-linker solution with $`\varphi `$ of order $`D/L`$ is expected to exhibit a sol-gel transition that is reasonably well described by percolation theory. The same argument indicates that the structure factor $`S(𝐪)`$ should resemble that of a pure Onsager nematic.

Recent high-precision low-angle synchrotron X-ray studies  of the microscopic properties of mixtures of actin filaments and positive divalent counterions report results that are in conflict with these expectations. When divalent counterions are added to a nematic phase of actin filaments, a birefringent gel phase is observed at linker concentrations that depends sensitively on the (average) rod length $`L`$, whereas percolation theory would predict that, beyond a certain length, the critical linker concentration is determined only by the rod concentration. In addition, a harmonic sequence is found in the structure factor $`S(𝐪)`$ of the gel phase in the long-wavelength range $`1/Lq1/D`$ along directions approximately perpendicular to the optical axis. This is incompatible with the Onsager nematic , but it is typical of layered structures, such as smectic liquid crystals. The X-ray experiments on actin gels as well as similar results obtained on other biopolymers such as DNA  clearly conflict with a simple “Onsager-Percolation” type description.

It is the claim of this paper, that the unusual features of biopolymer gels can be physically interpreted if we allow for the fact that linked biopolymers have preferential crossing angles. First, actin linker proteins are known to impose different crossing angles, depending on their molecular structure. Second, biopolymers like actin or DNA carry a high linear charge density, which allows them to be water soluble. The electrostatic repulsion between two uniformly charged rods linked at a single point depends on the relative angle $`\gamma `$ as $`1/|\mathrm{sin}\gamma |`$ , so the optimal crossing angle between two linked biopolymers is, in general, expected to be large. Based on this consideration, we studied \- as a simple model for biopolymer gel formation - a dilute mixture of rigid rods with a concentration $`\rho _p`$ of the order of $`1/DL^2`$, in the presence of a concentration $`\rho _l`$ of freely sliding, reversible linkers connecting pairs of rods at $`\pi /2`$ crossing angles (binding energy $`\epsilon _0`$). When $`\rho _l=0`$, second-virial theory predicts a first-order phase transition from an isotropic to a nematic phase. It is easy to demonstrate  that at non-zero linker concentrations, second-virial theory predicts a first-order phase transition (for $`\rho _l\rho _p`$) from the nematic phase to a birefringent “$`\pi /2`$” gel phase with tetragonal symmetry  (see Fig. 1b).

To examine gel formation beyond second-virial theory, we follow what happens when linkers are added, one by one, to the rod solution. A single linker connecting two rods produces a cross-like structure, with the two rods free to slide with respect to each other. It is easy to show that the free energy has a minimum when one arm of the cross extends along the optical axis (the $`z`$-direction) with the other arm taken to be the transverse $`x`$-direction . As shown in Fig. 2a, the transverse rod creates an anisotropic depletion volume $`\delta V_1`$ for the centers of mass (CM’s) of the free rods. This depletion volume depends on the angle $`\gamma _{ft}`$ between the axis of a free rod and the axis of the transverse rod by the Onsager result $`\delta V_1(\gamma _{ft})2DL^2|\mathrm{sin}\gamma _{ft}|`$. The entropic self-energy of the linked cross is computed as the osmotic work $`\delta W_1\mathrm{\Pi }_{\mathrm{osm}}\delta V_1(\gamma _{ft})_\theta `$ required to remove the CM’s of the free rods from the depletion volume, with $`\mathrm{\Pi }_{\mathrm{osm}}`$ the osmotic pressure of the nematic and with $`_\theta `$ denoting an orientational average. The effect of the arm of the cross along the optical axis is smaller by a factor $`D/L`$ and will be neglected. Using Onsager theory for the orientational distribution of the free rods, we find:

$$\frac{\delta W_1}{k_BT}\{\begin{array}{ccc}2c_I(1+c_I)& & (\mathrm{Isotropic},c_I3.3)\hfill \\ (24/\pi )c_N(1\mathrm{\Delta }\theta ^2/4)& & (\mathrm{Nematic},c_N4.5)\hfill \end{array}$$
(1)
with $`\mathrm{\Delta }\theta ^2`$ being the (small) variance of the polar angle with respect to the optical axis, and $`c=(\pi /4)DL^2\rho _p`$ the Onsager dimensionless rod concentration. At the isotropic/nematic transition point, the dimensionless free energy cost $`\delta W_1/k_BT33.2`$ assumes a large value independent of the dimensions of the rod, or the rod concentration.

Now add a second linker with a second rod making a $`\pi /2`$ angle with the optical axis (see example in Fig. 2b). Let the relative distance between the CM’s of the two transverse rods be $`𝐑`$, let the angle of the second rod with the $`x`$ direction (i.e., the relative angle) be $`\gamma _{tt}`$, and let $`\delta V_2(𝐑,\gamma _{tt})`$ be the excluded volume that is shared between the two rods. The entropic free energy cost of the two-rod system will then be

$$\delta W_2(𝐑,\gamma _{tt})2\delta W_1\mathrm{\Pi }_{\mathrm{osm}}\delta V_2(𝐑,\gamma _{tt})_\theta $$
(2)
In the limit $`R|𝐑|0`$ and $`\gamma _{tt}=0`$, the two rods coincide so $`\delta W_2(0,0)\delta W_1`$. On the other hand, for $`RL`$, there is no shared excluded volume so $`\delta W_2(𝐑,\gamma _{tt})2\delta W_1`$. It follows that the two linkers attract each other with a potential energy that has a typical range of order $`L`$ and a typical binding energy of order $`\delta W_1`$. This effect is related to the well-known depletion attraction , the attraction between large objects in a surrounding medium of small objects driven by the increase in configurational entropy as the larger objects coagulate. The difference is that in the present case an effective attraction between small objects (linkers) is generated by an increase in configurational entropy of large objects (rods). Using once again the Onsager orientational distribution, we find that, according to Eq. 2, the linker-linker depletion attraction drops off inversely proportional to the distance between the rods:

$$\frac{\delta W_2(z,\gamma _{tt})2\delta W_1}{\delta W_1}a_1\frac{D}{\mathrm{\Delta }\theta |z|}\mathrm{min}(\frac{D}{L|\mathrm{sin}\gamma _{tt}|},1)$$
(3)
where $`z`$ is the separation between the rods along the optical axis and $`a_16/\sqrt{\pi }`$. The interaction is highly anisotropic: for small $`|z|`$, there is a strong torque between transverse rods favoring parallel alignment (i.e., $`\gamma _{tt}=0`$). Moreover, the mechanism operates only as long as the CM of the second rod lies below or above the first rod within the two narrow wedges bordered by $`|y|\mathrm{\Delta }\theta |z|`$ centered on the first rod.

Equation 3 is valid for $`D\mathrm{\Delta }\theta |z|L`$. For small separations $`\mathrm{\Delta }\theta |z|D`$ the potential has a linear dependence on $`z`$:

$`{\displaystyle \frac{\delta W_2(z,0)2\delta W_1}{\delta W_1}}`$ (4)
$`\left(1a_2{\displaystyle \frac{\mathrm{\Delta }\theta |z|}{D}}\right)`$ $`\mathrm{min}({\displaystyle \frac{D}{L|\mathrm{sin}\gamma _{tt}|}},1)`$ (5)
with $`a_21/2\sqrt{\pi }`$. This result can be understood by noting that for perfectly aligned rods the shared depletion volume of two transverse rods would have a linear dependence on separation .

In the opposite limit of rod separations large compared to the rod length $`L`$, the continuum description of a nematic  can be used in terms of a locally varying director field $`\widehat{𝐧}(𝐫)`$ (see Fig. 3a). The perturbation of $`\widehat{𝐧}(𝐫)`$ far from the fixed rod can be written as $`\widehat{𝐧}(𝐫)(0,\delta n_y(𝐫),1)`$, where $`\delta n_y1`$. The leading term in the elastic free energy in this geometry is the splay term:

$$F_{\mathrm{splay}}=\frac{K_{11}}{2}\mathrm{d}^3r(\delta n_y)^2$$
(6)
where $`K_{11}`$ is the splay elastic constant and is equal to $`K_{11}0.06DL^4\rho _p^2k_BT`$ for rigid rods . Treating $`\delta n_y(𝐫)`$ as a variational function and minimizing Eq. 6 leads to the Laplace equation $`^2\delta n_y=0`$. As demonstrated in Fig. 3b, the appropriate solution of the Laplace equation must have the symmetry of a two-dimensional quadrupole so $`\delta n_y(y,z)=Qyz/(y^2+z^2)^2`$. The value of the quadrupole moment $`Q`$ is determined by the condition that the typical deviation $`Q/L^2`$ of $`\delta n_y(𝐫)`$ for $`y^2+z^2L^2`$ should match the angular deviation $`D/L`$ of the free rods as obtained from excluded volume arguments (see Fig. 2a). This condition leads to $`QDL`$ times a function of $`\mathrm{\Delta }\theta `$. From the electrostatic interaction energy between quadrupoles in two dimensions, $`Q^2\mathrm{cos}(4\phi )/r^4`$ (in polar coordinates $`r,\phi `$), it follows that the interaction energy per unit length between two parallel (infinite) rods must have the form:

$$\frac{F_{\mathrm{splay}}}{Lk_BT}g(\mathrm{\Delta }\theta )\mathrm{cos}(4\phi )c^2\frac{DL^2}{r^4}rL$$
(7)
where the function $`g(\mathrm{\Delta }\theta )`$ includes the dependence on the nematic order parameter. The deformation energy $`F_{\mathrm{splay}}`$ should be matched with the right hand side of Eq. 2 at distances $`rL`$ yielding $`g1/c`$. It follows from Eq. 7 that the anisotropy of the linker-linker interaction is significantly less pronounced in the regime $`rL`$. Moreover, the interaction is repulsive when the separation vector between the CM’s of the rods makes an angle near $`\pi /4`$ with the optical axis. The various regimes are shown schematically in Fig. 4.

We now turn to the implication of these results. According to Eq. 1, the characteristic energy scale $`\delta W_1`$ of the long-range depletion attraction between $`\pi /2`$ linkers is large compared to the thermal energy $`k_BT`$. This implies that percolation theory is not applicable since linkers will cluster even at very low concentrations. Next, second-virial theory is not reliable either since it gives no indication of this form of phase separation. Indeed, the cross-links resemble flat square discs of dimensions $`L\times L`$ and second-virial expansions are known not be accurate for mixtures of rods and discs . The large value of $`\delta W_1/k_BT`$ would have another important consequence: unless the linker binding energy is sufficiently high, crosses simply would not be able to form. Assuming chemical equilibrium between linkers in solution, linkers adsorbed on isolated chains, and linkers forming a cross-structure, we find that the fraction $`f`$ of linkers able to form a cross is of order:

$$f\frac{1}{1+(\pi /2c)\mathrm{e}^{(\epsilon _0\delta W_1)/k_BT}}$$
(8)
This means that cross-links disappear when $`\delta W_1`$ is large compared to $`\epsilon _0`$, i.e., when rods are so long that $`DL^2\epsilon _0/k_BT\rho _p`$. This is consistent with the experimental observation that the critical linker concentration for actin gelation rapidly increases with rod length. For an actin solution with $`\rho _p`$0.1mM, a linker concentration $`\rho _l`$1mM and $`\epsilon _015k_BT`$, the rod length should not exceed a value of order 100nm.

What would be the structure of the clusters? According to Eqs. 35, the depletion attraction between two transverse rods is highly anisotropic at short distances and maximal if the two rods are parallel and located just below each other with a common plane whose normal is perpendicular to the optical axis. Adding additional linkers and transverse rods in the same way produces an elongated ribbon structure resembling a disordered raft (see Fig. 1c). The sol-gel transition point should be characterized, not by percolation, but by the formation of elongated, ribbon-like rafts extending through the sample along the optical axis. Because two parallel plates immersed in a nematic attract each other  \- with a force per unit area that depends on the inter-layer spacing $`h`$ as $`k_BT/h^3`$ for large $`h`$ \- depletion attraction should also operate between different rafts, so the ribbons could be several layers thick (Fig. 1d). A multilayer ribbon-raft indeed would have along the transverse direction a structure factor $`S(𝐪)`$ similar to that of of a smectic, although without long-range order. Along the optical axis, there would have to be a peak in $`S(q_z)`$ when $`2\pi /q_z`$ equals the spacing between transverse rods inside the raft, which is indeed present .

Our method only applies to a small number of linkers. At finite linker concentration, rafts and ribbons are expected to act as long-lived kinetic intermediates. The actual thermodynamic equilibrium state is expected to be either a $`\pi /2`$ gel or a biaxial nematic , in phase-coexistence with linker-poor solution. Numerical simulation methods of the model should be able to address this issue. Next, we did not allow for the fact that linkers may prefer the endpoints of rods nor deviations away from $`\pi /2`$ for the preferential crossing angle. Finally, actin cytoskeletons do appear to resemble $`\pi /2`$ gels  with no evidence of raft formation. Actin filaments inside cells do not have a fixed length however: active polymerization and depolymerization processes are constantly taking place and in-vitro experiments of gel formation by such “living” actin filaments would be of great interest.

Acknowledgments: We would like to thank C. Safinya for communicating the results of the X-ray experiments on Actin gels as well as extensive discussions. We would like to thank P. van der Schoot and G.Wong for a critical reading of the manuscript and W. Gelbart, A. Liu, T. Lubensky and C. O’hern for useful discussions.


text:
# On fixed points of Poisson shot noise transforms

(Published in Advances in Applied Probability, 34, no.4, pp.798-825.)

## Abstract

Distributional fixed points of a Poisson shot noise transform (for nonnegative and nonincreasing response functions bounded by 1) are characterized. The tail behavior of fixed points is described. Typically they have either exponential moments or their tails are proportional to a power function, with exponent greater than $`1`$. The uniqueness of fixed points is also discussed. Finally, it is proved that in most cases fixed points are absolutely continuous, apart from the possible atom at zero.

Key words: Shot noise transform; fixed points; regular variation; renewal theorem; absolute continuity; infinite divisibility $``$ Banach Contraction Principle.
2000 Mathematics Subject Classification. Primary 60E07; Secondary 60K05


text:
# Unquenched Numerical Stochastic Perturbation Theory

## 1 Introduction

At Lattice 2000 we discussed how to include fermionic loops contributions in Numerical Stochastic Perturbation Theory for Lattice $`SU(3)`$, an algorithm which we will refer to as UNSPT (Unquenched NSPT). Our main message here is that unquenching NSPT results in not such a heavy computational overhead, provided only that an $`FFT`$ can be implemented in a fairly efficient way. $`FFT`$ is the main ingredient in constructing the fermion propagator by inverting the Dirac kernel order by order. For a discussion of the foundations of UNSPT we refer the reader to .

## 2 Lattice SU(3) UNSPT on APEmille

The need for an efficient $`FFT`$ is what forced us to wait for APEmille: our $`FFT`$ implementation mimic , which is based on a $`1dimFFT`$ plus transpositions, an operation which asks for local addressing on a parallel architecture. UNSPT has been implemented both in single and in double precision, the former being remarkably robust for applications like Wilson loops. To estimate the computational overhead of unquenching NSPT one can inspect Table 1. We report execution times of a fixed amount of sweeps both for quenched and unquenched NSPT. On both columns the growth of computational time is consistent with the the fact that every operation is performed order by order. On each row the growth due to unquenching is roughly consistent with a factor $`5/3`$. One then wants to understand the dependence on the volume, which is the critical one, the propagator being the inverse of a matrix: this is exactly the growth which has to be tamed by the $`FFT`$. One should compare execution times at a given order on $`L=8`$ and $`L=16`$ lattice sizes. Note that $`L=8`$ is simulated on an APEmille board ($`8`$ FPUs), while $`L=16`$ on an APEmille unit ($`32`$ FPUs). By taking this into account one easily understands that $`FFT`$ is doing its job: the simulation time goes as the volume also for UNSPT (a result which is trivial for quenched NSPT). Notice that at this level one has only compared crude execution times: a careful inspection of autocorrelations is anyway not going to jeopardize the picture. As for the dependence on $`N_f`$ (number of flavours), it is a parametric one: one plugs in various numbers and then proceed to fit the polynomial (in $`N_f`$) which is fixed by the order of the computation. It is then reassuring to find the quick response to a change in $`N_f`$ which one can inspect in Figure 1 (which is the signal for second order of the plaquette at a given value of the hopping parameter $`K`$).

## 3 Benchmark computation I: Wilson loops

We now proceed to discuss some benchmark computations. A typical one is given by Wilson loops. In Figure 2 one can inspect the first five orders 111All our expansions are written in powers of $`\beta ^1`$. of the basic plaquette at a given value of hopping parameter $`K`$, for which analytic results can be found in : going even higher in order would be trivial at this stage222Notice anyway that these computations are performed at a given value of hopping parameter $`K`$, but with no mass counterterm (see later).. Apart for being an easy benchmark, we are interested in Wilson loops for two reasons. First of all we are completing the unquenched computation of the Lattice Heavy Quark Effective Theory Residual Mass (see  for the quenched result). On top of that we also keep an eye on the issue of whether one can explain in term of renormalons the growth of the coefficients of the plaquette. There is a debate going on about that (see ), the other group involved having also started to make use of NSPT. In the renormalon framework the effect of $`N_f`$ can be easily inferred from the $`\beta `$-function, eventually resulting in turning the series to oscillating signs.

## 4 Benchmark computation II: the Wilson fermions Critical Mass

In Figure 3 we show the signal for one loop order of the Critical Mass for Wilson fermions (two loop results are available from ). The computation is performed in the way which is the most standard in Perturbation Theory, i.e. by inspecting the pole in the propagator at zero momentum. This is already a tough computation. It is a zero mode, an $`IR`$ mass-cutoff is needed and the volume extrapolation is not trivial. On top of that one should keep in mind that also gauge fixing is requested. The coefficients which are known analytically can be reproduced. Still one would like to change strategy in order to go to higher orders (which is a prerequisite of all other high order computations). The reason is clear: we have actually been measuring the propagator $`S`$, while the physical information is actually coded in $`\mathrm{\Gamma }_2=S^1`$ (one needs to invert the series and huge cancellations are on their way). Notice anyway that the fact that the Critical Mass is already known to two-loop makes many interesting computations already feasible.

## 5 Conclusions

Benchmark computations in UNSPT look promising, since the computational overhead of including fermionic loops contributions is not so huge. This is to be contrasted with the heavy computational effort requested for non perturbative unquenched lattice QCD. This in turn suggests the strategy of going back to perturbation theory for the (unquenched) computation of quantities like improvement coefficients and renormalisation constants. The Critical Mass being already known to two loops, many of these computations are already feasible at $`\alpha ^2`$ order.
We have only discussed the implementation of the algorithm on the APEmille architecture. We can also rely on a $`C^{++}`$ implementation for PC’s (clusters) which is now at the final stage of development.


text:
# Quintessence, Cosmology, and FRIIb Radio Galaxies

## 1 Introduction

There are several independent ways to determine the global cosmological parameters that describe the current state of the universe. It is important to have several complementary and independent methods that yield consistent results since any given method could be plagued by unknown systematic errors. A particularly useful way to determine global cosmological parameters is through measurements of the coordinate distance to high-redshift sources. This method is particularly useful because the coordinate distance depends only upon global cosmological parameters, and is independent of the clustering properties of the mass-energy components that control the expansion rate of the universe, as long as each component is homogeneous and isotropic on very large scales. The coordinate distance is also independent of whether different components cluster differently, known as biasing, and on whether the dark matter is cold, warm, or hot (though it does depend on the equation of state of each component).

Two cosmological tools that are particularly sensitive to the coordinate distance as a function of redshift are FRIIb radio galaxies, which provide a modified standard yardstick (Daly 1994, 2002a; Guerra 1997; Guerra & Daly 1998 \[GD98\]; Guerra, Daly, & Wan 2000 \[GDW00\]), and type Ia supernovae, which provide a modified standard candle (e.g. Riess et al. 1998, Perlmutter et al. 1999). Some aspects of the methods are compared in §4.

The coordinate distance to a source at a given redshift depends on the present value of the mean mass-energy density of each component, and on the redshift evolution of the mass-energy density of each component. Non-relativistic matter, including baryonic matter and clustered dark matter, have a mean mass-energy density that evolves as $`(1+z)^3`$. There is a contribution from radiation and neutrinos left over from the big bang, currently negligible, which has an energy-density that evolves as $`(1+z)^4`$. There may also be a cosmological constant, which constant energy density, and so evolves as $`(1+z)^0`$.

The redshift evolution of the mean mass-energy density of a component depends on the equation of state $`w=P/\rho `$ of the component. There could exist a dynamical vaccum energy density called Quintessence (Caldwell, Dave, & Steinhardt 1998). Quintessence would have an equation of state $`w`$ such that the density evolves as $`(1+z)^n`$, where $`n=3(w+1)`$ (see Turner & White 1997; Bludman & Roos 2001, Wang et al. 2000; or the Appendix of this paper). For example, pressure-less dust (i.e., non-relativistic matter) is described by $`P=w=0`$ and $`n=3`$; a relativistic fluid is described by $`w=1/3`$ and $`n=4`$, and a cosmological constant is described by $`w=1`$ and $`n=0`$.

Two components of the universe are known to exist with certainty: there is a non-relativistic component (including baryons and clustered dark matter), with zero-redshift, normalized mean mass-energy density $`\mathrm{\Omega }_m`$, and a primordial relativistic component that is negligible for the epochs of interest here (including microwave background photons and neutrinos).

There must be at least one more component, which remains to be determined. This third component could be a cosmological constant, quintessence, space curvature, or something else. In the simplest case, there is only one additional component that is significant and has a measurable impact over the redshift interval where the coordinate distance is used to constrain cosmological parameters and the properties of the unknown component.

Here, the determination of the coordinate distances to 20 radio galaxies with redshifts between zero and two are used to constrain the equation of state and current mean mass-energy density of quintessence assuming a spatially flat universe as indicated by measurements of the microwave background radiation (de Bernardis et al. 2000, Balbi et al. 2000).

The radio galaxy method is briefly reviewed in §2. Constraints on cosmological parameters and the equation of state of quintessence obtained using FRIIb radio galaxies as a modified standard yardstick are presented in §3. The supernovae and radio galaxy methods and results are compared in detail in §4. The relation between beam power and total energy for FRIIb radio jets is discussed in §5. Here it is shown that the underlying hypotheses of the radio galaxy method is consistent with current models of large-scale jet production in AGN. The results are summarized in §6.

## 2 Cosmology FRIIb Radio Galaxies

The use of FRIIb radio galaxies as a modified standard yardstic to determine global cosmological parameters has been presented in detail elsewhere; see Daly (1994, 2002a), GD98, and GDW00. Radio images of the 20 sources used here can be found in Leahy, Muxlow, & Stephens (1989), Liu, Pooley, & Riley (1992), and GDW00. The radio galaxy method is based on the following premises. First, it is assumed that the forward region of the radio bridge of an FRIIb source represents a strong shock front, and the equations of strong shock physics are applied at this boundary. Second, it is assumed that there is a power-law relation between the average beam power of the source $`L_j`$ and the total lifetime of the outflow $`t_{}`$: $`t_{}L_j^{\beta /3}`$ (see §5). Third, it is assumed that, at a given redshift, the sources studied will have an average source size similar to that of the parent population of FRIIb sources at the same redshift, where the source size is defined as the separation between the radio hot spots. All of these points are discussed in detail in previous work (e.g. Daly 2002a,b; GDW00).

In addition, in its simplest interpretation, the model works in a straight-forward way if the radio power of a given FRIIb radio galaxy does not decrease monotonically and measurably with time. There is substantial evidence that this is the case for radio galaxies that fall into the FRIIb category, as do all of the radio galaxies considered here (Daly 2002b). For example, Neeser et al. (1995) find that there is no correlation between radio size and radio power. The low frequency radio surface brightness of these sources can be explained using adiabatic expansion alone if the radio powers of the sources are roughly time-independent (Wellman, Daly, & Wan 1997). Gopal-Krishna, Kulkarni, & Wiita (1996) find that a roughly constant radio power followed by an exponential drop in radio power would explain the relative sizes of radio galaxies and radio loud quasars in the context of the orientation unified model. And, it is shown in §5 that the beam power is independent of core-hot spot separation, which would not be expected if the radio power of a given source decreases with time. Thus, there is no empirical evidence that the radio power of a given FRIIb radio galaxy decreases with time, and substantial evidence that it is roughly constant over the lifetime of a given source, followed by a precipitous drop. The Mach numbers of FRIIb radio galaxies are about 2.5 to 10 (Wellman, Daly, & Wan 1997), implying overpressures of 10 to 200, so the lobe will rapidly expand and drop in radio power when the beam power shuts off, as suggested by Neeser et al. (1995). Some authors (e.g., Blundell, Rawlings, & Willcot 1999) have studied large heterogeneous samples that include sources with all types of radio structure and power and both radio loud quasars and radio galaxies. To explain the properties of this heterogenous sample, a model is proposed that allows for the decrease in the radio power of a given source as a function of time. There is, however, no evidence that this applies to the select subset of sources (the most powerful FRII radio galaxies) studied here, and substantial evidence that the radio powers of FRIIb radio galaxies do not monotonically decline over the lifetime of a given source.

Since FRIIb radio galaxies form a remarkably homogenous population, and all the sources at a given redshift are quite similar, the average size of a given source $`D_{}`$ will equal the average size of the full population $`<D>`$ of FRIIb radio galaxies with similar redshifts. The dependence of $`D_{}`$ and $`<D>`$ on the coordinate distance, and hence on cosmological parameters, is quite different, and

$$<D>/D_{}(a_or)^{2\beta /3+3/7}$$
(1)
since $`<D>(a_or)`$ (see GD98 eq. 6). By construction, the ratio $`<D>/D_{}`$ should be a constant, independent of redshift. Because the model is based upon a comparison of individual source characteristics with the characteristics of the parent population at the same redshift, selection effects, such as radio power selection effects, are unlikely to affect the constraints placed on cosmological parameters. All of the sources for which values of $`D_{}`$ are determined belong to the parent population. The sources for which $`D_{}`$ is determined are subject to selection effects that are identical to those of the parent population, used to determine $`<D>`$. Thus, this method is one of the few methods that is unlikely to be affected by observational selection effects. Equation (1) allows constraints to be placed on the model parameter $`\beta `$, and the cosmological parameters that enter into the determination of the coordinate distance.

## 3 Constraints on Quintessence From Radio Galaxies

Recent measurements of the cosmic microwave background anisotropy suggest that the universe has zero space curvature (de Bernardis et al. 2000, Balbi et al. 2000). Here, a universe with zero space curvature ($`k=0`$) and quintessence is considered. The coordinate distance to a source at redshift $`z`$ in a spatially flat universe is

$$(a_or)=H_o^1_0^z𝑑z/E(z).$$
(2)
where $`E(z)=\sqrt{\mathrm{\Omega }_m(1+z)^3+(1\mathrm{\Omega }_m)(1+z)^n}`$. The coordinate distance is obtained by substituting this into eq. 16 and solving for $`(a_or)`$ (see, for example, Turner & White 1997; Bludman & Roos 2001; the Appendix of this paper).

The equation of state $`w`$ of quintessence is assumed to be time-independent over the redshift interval from zero to two (Wang et al. 2000). Thus, the mean mass-energy density of quintessence will evolve as $`(1+z)^n`$, where $`w=n/31`$. For a present value of the normalized mass-energy density of quintessence $`\mathrm{\Omega }_Q`$, and that of non-relativistic matter $`\mathrm{\Omega }_m`$, we have $`\mathrm{\Omega }_Q+\mathrm{\Omega }_m=1`$ at $`z=0`$, since space curvature is taken to be zero. This implicitly assumes that there is only one unidentified component, $`\mathrm{\Omega }_Q`$, that contributes significantly over the redshift interval where the coordinate distance will be determined.

The minimization of the chi-squared of $`<D>/D_{}=`$ constant leads to constraints on $`\mathrm{\Omega }_m`$, $`w`$,and the model parameter $`\beta `$. The one and two-dimensional confidence intervals have been determined in the usual way (see, for example, Press, Teukolsky, & Flannery 1992).

Figure 1 shows the 68%, 90%, and 95% one-dimensional confidence contours in the $`\mathrm{\Omega }_mw`$ plane. At 90% confidence, $`\mathrm{\Omega }_m`$ is less than about 0.5 for any value of $`w`$. A cosmological constant, represented by $`w=1`$ is clearly consistent with the radio data. This Figure differs slightly from the two-dimensional confidence intervals presented in Figure 2 of Daly & Guerra 2001. At 68 % confidence, $`\mathrm{\Omega }_m<0.25`$, and $`1.3<w<0.45`$. At 90 % confidence, $`\mathrm{\Omega }_m<0.5`$ and $`2.6<w<0.25`$. Note, some models for the dark energy invoke mechanisms that behave like quintessence with w less than -1 (e.g.Frampton 2002).

The fifty-four supernovae identified by Perlmutter et al. (1999) and used to obtain their “primary fit” C have been similarly fitted, and imply strong contraints on $`\mathrm{\Omega }_m`$, though the supernovae do not place a lower bound on w, as shown in Figure 2.

As described by Turner & White (1997) and presented in some detail in the Appendix, the expansion rate of the universe is accelerating at the present epoch if

$$1+3w(1\mathrm{\Omega }_m)<0.$$
(3)
This line is drawn on the $`\mathrm{\Omega }_mw`$ plane (see Figures 1 and 2); points below the line represent solutions for which the universe is currently accelerating, while those above the line represent solutions for which the universe is currently decelerating. Radio galaxies alone (Fig 1) indicate that the expansion of the universe is accelerating ($`q_o<0`$); this is an 84% confidence result.

It is very important to consider whether significant covariance exists between the different parameters determined by the fit. Figure 3 shows the 68%, 95%, and 99% one-dimensional confidence contours in the $`w\beta `$ plane. Clearly, there is no covariance between the equation of state $`w`$ of quintessence and the model parameter $`\beta `$: $`w`$ and $`\beta `$ are independent.

Figure 4 shows the 68%, 95%, and 99% one-dimensional confidence contours in the $`\mathrm{\Omega }_m\beta `$ plane. Clearly the two parameters are independent. The data suggest that $`\beta 1.75\pm 0.25`$ independent of the cosmological parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_Q`$, or the equation of state $`w`$ of quintessence.

## 4 Comparison of Supernova and Radio Galaxy Methods

The observation of the acceleration of the universe through measurements of the coordinate distance to high redshift sources is so important that it needs to be independently verified using different methods. In §4, it was shown that the radio galaxy method alone indicates at about 84% confidence that the universe is accelerating in its expansion at the present epoch if the universe has zero space curvature. The supernovae teams report similar results at stronger confidence levels (Riess et al. 1998, Perlmutter et al. 1999), as illustrated in Figure 2. It is important to consider the similarities and differences between the methods to insure that the results of each are independent, and to compare the quantities that are actually determined in each method. The two methods are compared in Table 1.

The two methods are similar is their dependence on the coordinate distance, and in the sense that both provide a modified standard rather than an absolute standard. They are somewhat similar in their redshift interval of coverage, and in the number of sources studied. The radio galaxy and supernova methods differ in the way the “modification” relation is derived; one is empirical and one is based on physical arguments. They also differ in that the radio galaxy method is not normalized at zero redshift, and is independent of the properties of local sources and of the local distance scale; in fact, nearly identical results obtain when the low redshift bin is excluded from the analysis (GDW00). The supernova relation is derived at low redshift.

Thus, the radio galaxy method provides independent verification of the supernova results. Current constraints on cosmological parameters determined using the radio galaxy method are consistent with those obtained using the supernova method. Clearly, any selection effects, or unknown systematic errors, must be completely different for these two methods. The supernova method involves short-lived optical events, and the radio galaxy method involves much more long-lived radio AGN events. The methods are derived in vastly different ways, and are normalized differently. However, similar constraints obtain, suggesting each method is reliable.

### 4.1 Comparison Allowing for Space Curvature

It is interesting to compare the constraints placed on cosmological parameters assuming that three main terms control the expansion of the universe at the present epoch: non-relativistic matter with contribution $`\mathrm{\Omega }_m`$, a cosmological constant with contribution $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and space curvature. Figures 5 and 6 show constraints placed by the microwave background radiation (Bond et al. 2000), and the high-redshift supernovae project (Riess et al. 1998); the results obtained by the supernovae cosmology project (Perlmutter et al. 1999) are similar and are presented by Daly & Guerra (2001). Each method carves out a different part of parameter space because each method covers a different redshift range with the supernovae being at the lowest redshift, the radio galaxies at higher redshift, and the microwave background radiation at still higher redshift. The three methods are thus complementary in their redshift coverage.

The radio galaxy and supernovae methods are independent of the initial power spectrum of density fluctuations, and hence are independent of the index of the primordial fluctuation spectrum. These methods are also independent of the Hubble constant, and of the baryon fraction. The radio galaxy and supernovae methods depend only on the global cosmological parameters $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ allowing for space curvature, or $`\mathrm{\Omega }_m`$ and $`w`$ allowing for quintessence, where $`\mathrm{\Omega }_Q=1\mathrm{\Omega }_m`$, and assuming zero space curvature. They are therefore particularly important to study since most other methods depend on several other parameters such as the initial spectrum of fluctuations, the Hubble constant, and/or the baryon fraction.

## 5 Lifetime, Beam Power, and Energy in FRIIb Radio Galaxies

The application of the radio galaxy cosmology method indicates that $`\beta 1.75\pm 0.25`$ irrespective of cosmological model, as discussed here for quintessence, by GDW00 in cosmologies that include space curvature, and summarized by Daly (2002a). The total energy released over the lifetime of the FRIIb in the form of a highly collimated outflow is $`E_{}=L_jt_{}`$ assuming that $`L_j`$ is roughly constant over the lifetime of a given source, which is supported by lack of correlation between $`L_j`$ and the hotspot to core distances in these sources (see Wan, Daly, & Guerra 2000; Guerra & Daly 2001). Values of $`L_j`$ and $`t_{}`$ are obtained independently, using different aspects of the radio data for a given source. Two values are determined for each source, one on each side of each source. An indicator of possible problems with the model assumptions or selection biases is whether the beam powers, total source lifetimes, or total energies depend strongly on the distance of the radio hot spot from the location of the AGN. Each of these quantities is plotted as a function of hot spot distance from the AGN in Figures 7, 8, and 9. There are 40 data points on each figure, two for each of the 20 radio galaxies studied here.

It is interesting to note that the total energy processed through large-scale jets has a rather small range, from $`few\times 10^5M_{}`$ to $`few\times 10^6M_{}`$, though it does increase systematically with redshift, as shown in Figure 10.

Because $`t_{}L_j^{\beta /3}`$ and $`E_{}=L_jt_{}`$, the relation between beam power and energy released by an FRIIb is $`L_jE_{}^{3/(3\beta )}E_{}^q`$, where $`q=2`$ to 3 for $`\beta =1.5`$ to 2.0. In terms of the total lifetime of the source, $`t_{}L_j^{(1/2,7/12,2/3)}`$ for values of $`\beta `$ of 1.5, 1.75, and 2.0 respectively. Clearly, $`\beta `$ remains far from the limiting cases of $`\beta `$ = 0 or 3.

### 5.1 Electromagnetic Energy Extraction from Rotating Holes

Blandford (1990) gives a summary of the energy and beam power of outflows associated with the electromagnetic extraction of the spin energy of a rotating black hole. The beam power and total energy available are

$$L_j=L_{EM}10^{45}(a/m)^2B_4^2M_8^2\text{ erg s}^1(a/m)^2B^2M^2$$
(4)
and

$$E_{}=E5\times 10^{61}(a/m)^2M_8\text{ erg}(a/m)^2M,$$
(5)
for $`(a/m)<<1`$, where $`M`$ is the mass of the black hole, $`M_8`$ is the mass in units of $`10^8M_{}`$, $`a`$ is the spin angular momentum $`S`$ per unit mass $`M`$: $`a=S/(Mc)`$, c is the speed of light, $`m`$ is the gravitational radius $`m=GM/c^2`$, $`B`$ is the magnetic field strength, and $`B_4`$ is the magnetic field strength in units of $`10^4`$ G (Blandford 1990).

Empirically, $`E_{}L_j^{1\beta /3}`$ with $`\beta =1.75\pm 0.25`$. This is consistent with equations (4) and (5) when the magnetic field strength satisfies

$$BM^{(2\beta 3)/2(3\beta )}(a/m)^{\beta /(3\beta )}.$$
(6)
This is particularly simple when $`\beta =1.5`$; in this case $`B(a/m)`$. For the cases $`\beta =1.75`$ and 2, the empirical constraint is consistent with equations (4) and (5) when $`B(a/m)^{1.4}M^{0.2}`$ and $`B(a/m)^2M^{1/2}`$ respectively.

Typical values for the beam power, total energy processed by a jet, and lifetime are $`L_j10^{45}\text{ erg s}^1`$, $`E_{}5\times 10^5M_{}c^2`$, and $`t_{}10^7`$ yr ( see Figures 7, 8, and 9). Equations 4 and 5 imply a total source lifetime of $`t_{}10^9/(M_8B_4^2)`$, indicating a lifetime of about $`10^7`$ yr for $`M_8B_4^210^2`$. This is satisfied for $`M_810`$, and $`B_43`$. In this case, the beam power is $`10^{45}\text{ erg s}^1`$ for $`(a/m)(1/30)`$, and the total energy is $`5\times 10^5M_{}`$. These values for $`M_8`$, $`B_4`$, and $`(a/m)`$ seem quite reasonable. The scaling between variables required by the empirically determine relation between total energy and beam power are discussed above.

## 6 Conclusion

Type IIb radio galaxies and type Ia supernovae are particularly important methods to develop to constrain the global cosmological parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and space curvature, or $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_Q`$, and the equation of state of quintessence $`w`$. These methods depend only upon the properties of global cosmological parameters, and are independent of other factors such as the index of the primordial power spectrum, the Hubble constant, the baryon fraction, the properties of the dark matter that clusters around galaxies and clusters of galaxies, any biasing of dark relative to luminous matter, etc. The radio galaxy and supernova methods are completely independent and have completely different potential systematic errors, as discussed in §4.

Radio galaxies may be used to constrain the mass-energy density and equation of state of quintessence. These results are presented here assuming a spatially flat universe, which is supported by recent measurements of fluctuations of the cosmic microwave background. Radio galaxies alone suggest that the universe is accelerating in its expansion at present (see S 3), consistent with results obtained by the supernovae teams, indicating that both methods are working well and probably are not plagued by unknown systematic errors.

The implications for models of energy extraction for the cases of an outflow related to accretion and the Eddington luminosity, and electromagnetic energy extraction of rotational energy are considered. If the outflow is produced by the electromagnetic extraction of energy from a rotating black hole, then the magnetic field strength must be related to the spin angular momentum of the rotating black hole $`S`$, the mass of the black hole $`M`$, and the gravitational radius of the black hole $`m`$, as described in §5. The relation is particularly simple if $`\beta =1.5`$, and implies that the magnetic field strength satisfies $`B(a/m)`$, where $`a=S/(Mc)`$.

We are grateful to Lin Wan and Greg Wellman who have contributed significantly our studies of radio galaxies. We also thank Megan Donahue, Paddy Leahy, Chris O’Dea, Adam Riess, Max Tegmark, and an anonymous referee for helpful comments. This work was supported in part by a National Young Investigator Award AST-0096077 from the National Science Foundation, and by Penn State University. Research at Rowan University was supported in part by the College of Liberal Arts and Sciences and National Science Foundation grant AST-9905652.

## 7 Appendix

The deceleration parameter is defined to be $`q_o=(\ddot{a}a)/\dot{a}^2`$ evaluated at z=0, where $`a`$ is the cosmic scale factor. This may also be written $`q_o=\ddot{a}_o/(a_oH_o^2)`$, where quantities evaluated at z=0 have a subscript ’o’, and $`H_o=(\dot{a}_o/a_o)`$. It is straight-forward to show that (e.g. Peebles  eq. 5.15)

$$(\ddot{a}/a)=\frac{4\pi G}{3}(\rho _i+3p_i)=\frac{4\pi G}{3}\rho _i(1+3w_i)$$
(7)
where $`p_i`$ is the pressure, $`\rho _i`$ is the mean mass-energy density, $`w_i`$ is the equation of state of $`i\text{th}`$ component, $`w_i=p_i/\rho _i`$, and the quintessence or cosmological constant term is included in the summation.

For each component, mass-energy conservation implies

$$\dot{\rho _i}=3(\rho _i+p_i)(\dot{a}/a)$$
(8)
(e.g. Peebles  eqs. 5.16). By definition, the equation of state is $`w_i=p_i/\rho _i`$, so the eq. 8 implies $`(\dot{\rho _i}/\rho _i)=3(1+w_i)(\dot{a}/a).`$ When the equation of state $`w_i`$ is time independent, the solution to this equation is $`\rho _i=\rho _{i,o}(1+z)^{3(1+w_i)}`$, where $`(1+z)=a_o/a`$. Thus, a component with equation of state $`w_i`$ and present mean mass-energy density $`\rho _o`$ will have a mean mass-energy density at redshift $`z`$ of $`\rho =\rho _o(1+z)^{n_i}`$, where $`n_i=3(1+w_i)`$.

When $`k=0`$,

$$(\dot{a}/a)^2=\frac{8\pi G}{3}\rho _i,$$
(9)
where $`\rho _i=\rho _{i,o}(1+z)^n`$. It follows that $`H_o^2=(\dot{a}_o/a_o)^2=(\frac{8\pi G}{3})\rho _{i,o}`$. For $`k=0`$, $`\rho _{i,o}=\rho _{c,o}`$ where $`\rho _{c,o}`$ is the critical density at redshift zero. Thus,

$$H_o^2=\frac{8\pi G}{3}\rho _{c,o}$$
(10)
for $`k=0`$. The deceleration parameter $`q_o=\ddot{a}_o/(a_oH_o^2)`$ then becomes

$$q_o=(1/2)\mathrm{\Omega }_i(1+3w_i),$$
(11)
where $`\mathrm{\Omega }_i=\rho _{i,o}/\rho _{c,o}`$ and equations 7 and 10 have been used.

The universe is accelerating rather than decelerating when $`q_o<0`$. This can only occur if $`(1+3w_i)<0`$, or $`w_i<1/3`$, which is a necessary but not sufficient condition to have an accelerating universe at the present epoch. If there are only two significant types of mass-energy controlling the expansion rate of the universe at the present epoch, quintessence and non-relativistic matter, then the deceleration parameter is $`q_o=\mathrm{\Omega }_m/2+\mathrm{\Omega }_Q(1+3w)/2`$. The universe will be accelerating in its expansion when

$$1+3w(1\mathrm{\Omega }_m)<0,$$
(12)
which follows since $`\mathrm{\Omega }_Q=1\mathrm{\Omega }_m`$. Thus, a curve can be drawn on the $`\mathrm{\Omega }_mw`$ plane; points below the curve represent solutions for which the universe is currently accelerating, while those above the curve represent solutions for which the universe is currently decelerating (e.g. see Figure 1). For example, this curve is shown on Figure 1, which illustrates the constraints obtained using radio galaxies alone, and is discussed in more detail below. Note that radio galaxies alone place interesting constraints on $`\mathrm{\Omega }_m`$ and $`w`$, and these results are consistent with those obtained using other methods (e.g Wang et al. 2000).

The coordinate distance to a source at redshift $`z`$ follows from the equation

$$𝑑r/\sqrt{1kr^2}=𝑑t/a(t)=(1/a_o)(\dot{a}/a)^1𝑑z$$
(13)
(Weinberg 1972). For a spatially flat universe, the left hand side of the equation reduces to the coordinate distance $`r`$. Equations 9 and 10 imply that $`(\dot{a}/a)^2=(H_o^2/\rho _{c,o})\rho _{i,o}(1+z)^{n_i}`$ or

$$(\dot{a}/a)^2=H_o^2\mathrm{\Omega }_i(1+z)^{n_i}.$$
(14)
Following the notation of Peebles (1993), this reads $`(\dot{a}/a)=H_oE(z)`$, where

$$E(z)=\sqrt{\mathrm{\Omega }_i(1+z)^{n_i}}.$$
(15)
Thus, the coordinate distance to a source at redshift $`z`$ in a spatially flat universe is (see equation 13)

$$(a_or)=H_o^1_0^z𝑑z/E(z).$$
(16)
A spatially flat universe has $`\mathrm{\Omega }_i=1`$. The components that must be included in equation 15 are those that contribute from redshift zero out to the redshift that the coordinate distance is being determined. Two components are considered here: non-relativistic matter with normalized mean mass-energy density at $`z=0`$ of $`\mathrm{\Omega }_m`$, and quintessence with normalized mean mass-energy density at $`z=0`$ of $`\mathrm{\Omega }_Q`$. Thus, $`\mathrm{\Omega }_m+\mathrm{\Omega }_Q=1`$, and equation 15 becomes $`E(z)=\sqrt{\mathrm{\Omega }_m(1+z)^3+(1\mathrm{\Omega }_m)(1+z)^n}`$. The coordinate distance is obtained by substituting this into eq. 16 and solving for $`(a_or)`$.


text:
# Defect correlations, metal-insulator transition, and magnetic order in ferromagnetic semiconductors

## 1 Introduction

Diluted ferromagnetic III-V semiconductors are promising materials for applications as well as interesting from the physics point of view. They might allow the incorporation of ferromagnetic elements into semiconductor devices, and thus the integration of data processing and magnetic storage on a single chip. It is important to understand the interplay between transport properties, magnetic ordering, and the defect configuration, which is determined by the growth conditions. We show that correlated positions of defects are required for a description consistent with experiments .

To be specific, we consider GaAs doped with manganese, but the results should apply for most III-V compounds. The physics relies on the dual role played by the Mn impurities: They carry a local spin due to the half-filled *d*-shell and dope the system with holes, which mediate a ferromagnetic indirect exchange interaction between the spins. A crucial feature of these materials is their high degree of compensation, which is presumably due to antisite defects (As substituted for Ga). The accordingly small density of holes leads to weak electronic screening of the Coulomb interaction between charged defects and between defects and holes. This role of compensation has not been taken into account previously. Due to the strong Coulomb interactions a random distribution of defects on the cation sublattice is very costly in energy. Therefore, defect diffusion  leads to a rearrangement of defects in such a way that the Coulomb energy is reduced. We study this rearrangement with the help of Monte Carlo (MC) simulations and find that strong correlations of defect positions develop. The resulting disorder potential acting on the valence-band holes is strongly reduced and most of its correlations fall off on the scale of the minimum defect separation. We show that this effect has tremendous consequences for transport properties, the integrity of the energy gap, and the temperature-dependent magnetisation. In fact, a description in agreement with experiments *requires* the defects to be strongly correlated.

## 2 Correlated defects

In this section, we study how diffusion of defects, *i.e.*, Mn impurities and antisites, changes the defect configuration and the disorder potential. This is motivated by experiments  that show that defect diffusion is rather pronounced at typical growth and annealing temperatures of the order of $`250^{}`$C. We start from a random distribution of Mn impuritites and antisites (As) on the cation sublattice, where the density of antisites is determined by charge neutrality from the observed  hole concentrations at given Mn doping levels. The Hamiltonian reads

$$H=\frac{1}{2}\underset{i,j}{}\frac{q_iq_j}{ϵr_{ij}}e^{r_{ij}/r_{\mathrm{scr}}},$$
(1)
where $`q_i`$ are the defect charges and $`r_{ij}`$ is their separation. Relative to the cation sublattice, Mn impurities carry charge $`q=1`$ whereas antisites carry $`q=+2`$. The screening length $`r_{\mathrm{scr}}`$ is obtained from nonlinear screening theory . It is much larger than the nearest-neighbor separation on the cation sublattice so that it hardly affects the small-scale defect correlations relevant here. We perform MC simulations for the Hamiltonian $`H`$  employing the Metropolis algorithm at the temperature $`250^{}`$C for systems of $`20\times 20\times 20`$ conventional cubic unit cells with periodic boundary conditions, unless stated otherwise. To roughly model the dynamics, we only use local MC moves, *i.e.*, exchanges of nearest-neighbor defect and/or host atoms. A more realistic modelling of the dynamics involving interstitials and vacancies lies beyond the scope of this work. The resulting configurations are quenched at low temperatures .

Figure 1 (left) shows a random configuration of defects for a Mn concentration of $`x=0.05`$ and $`p=0.3`$ holes per Mn, which determines the concentration of antisites ($`0.0175`$). This should be compared to the right part of Fig. 1, which shows an equilibrated configuration of the same number of defects. It is obvious that the defects have formed small clusters, where typically an antisites (charge $`+2`$) is surrounded by several Mn impurities (charge $`1`$). It is clear that this clustering leads to a strongly reduced impurity potential, since the antisite charges are very effectively screened.

More quantitative information can be obtained from the potential correlation function

$$D(r)V(𝐫)V(𝐫^{})_{\left|𝐫𝐫^{}\right|=r}V^2.$$
(2)
Obviously, $`\mathrm{\Delta }V\sqrt{D(0)}`$ is the width of the distribution of $`V(𝐫)`$. $`D(r)`$ is plotted in Fig. 2 for $`x=0.05`$ and $`p=0.3`$ for various numbers of MC steps, which serve as a rough meassure of annealing time. Defect diffusion strongly reduces $`D(r)`$ and thus $`\mathrm{\Delta }V`$. The inset in Fig. 2 shows a blowup of $`D(r)`$ for the approximately equilibrated state. The initial reduction of $`D(r)`$ to about 10% of $`D(0)`$ on the scale of the nearest-neighbor separation of about $`4.00`$Å(heavy bar) is due to the screening of the compensated 70% of Mn impurities by antisites, leaving only 30% active. Since $`D(r)`$ contains the potential squared, this leads to a reduction to about 9% (yellow region in the inset). The remaining contribution from uncompensated Mn (light blue region) decays on the typical length scale of $`14.4`$Åassociated with their density, which is $`0.3\times 0.05`$ per cation site. Clearly, ionic screening, *i.e.*, screening by the defects, is nearly perfect at $`250^{}`$C. While the reduction of $`\mathrm{\Delta }V`$ is thus substantial, we find that $`\mathrm{\Delta }V`$ is still *not* small compared to the Fermi energy .

## 3 Valence-band holes

In this section we consider the properties of valence-band holes  in the impurity potential $`V(𝐫)`$ due to the charged defects. In particular, we are interested in their spectrum and localisation properties. We employ the envelope function and parabolic-band approximations for the holes and start from the Hamiltonian $`H=_i(\mathrm{}^2/2m^{})_i^2+V(𝐫_i)`$. For material-specific calculations the detailed band structure should be taken into account, *e.g.*, using the 6-band Kohn-Luttinger Hamiltonian . This should not change the qualitative results but is known to increase the mean-field $`T_c`$ . The hole Hamiltonian is written in a plane-wave basis and diagonalized numerically . The calculations are done for spin-less holes, since the additional disorder introduced by the exchange interaction is found to be much smaller than $`\mathrm{\Delta }V`$. We obtain the energy spectrum and normalized eigenfunctions $`\psi _n(𝐫)`$. From the latter we calculate the inverse participation ratios

$$\mathrm{IPR}(n)=1/\underset{𝐫}{}|\psi _n(𝐫)|^4$$
(3)
of the states. The main physical content of the IPR is that it scales with system size for extended states but essentially remains constant for localized states. The IPR thus allows to estimate the position of the mobility edge in the valence band.

Figure 3 shows the IPR as a function of electron energy for the same parameters as above, $`x=0.05`$ and $`p=0.3`$, after various numbers of MC steps. The plot shows that uncorrelated defects would lead to the filling-in of the band gap by disorder, which is in contradiction to experiments. Thus correlated defects are required to explain the persistence of the energy gap. We have also studied the IPR for various system sizes and find that states in the flat region and the upper part of the slope are extended, whereas states in the band tail are localized, as expected. Figure 3 shows that the states at the Fermi energy become more and more extended with annealing, since the disorder potential decreases.

In Fig. 4 the IPR for fully annealed configurations is plotted as a function of Mn concentration $`x`$. The number of holes per Mn, $`p`$, has been chosen in accordance with experiments . The change in the IPR-curves with $`x`$ is rather weak due to the strong ionic screening. The main effect is the shift of the Fermi energy due to the varying hole density, which drives the Fermi edge into localized states below about 1% Mn. This qualitatively explains the metal-insulator transition observed in (Ga,Mn)As .

## 4 Ferromagnetism

Finally, we turn to the magnetic properties. The bulk magnetisation is composed of contributions of the local Mn spins and of the holes. Since the exchange interaction between valence-band holes and local Mn spins is antiferromagnetic, these contributions are oriented oppositely. However, the hole contribution turns out to be very small. We study the magnetisation within a selfconsitent mean-field approximation, which we briefly outline in the following.

The starting point is the Hamiltonian of the coupled valence-band holes and Mn spins,

$$H=\underset{n\sigma }{}\xi _nc_{n\sigma }^{}c_{n\sigma }J_{\mathrm{pd}}\underset{i}{}𝐬_i𝐒_i,$$
(4)
where $`𝐒_i`$ are Mn spins ($`S=5/2`$) and

$$𝐬_i\underset{n\sigma n^{}\sigma ^{}}{}c_{n\sigma }^{}\psi _n^{}(𝐑_i)\frac{𝝈_{\sigma \sigma ^{}}}{2}\psi _n^{}(𝐑_i)c_{n^{}\sigma ^{}}$$
(5)
are hole spin polarisations at the Mn sites. $`\xi _n`$ and $`\psi _n`$ are the hole eigenenergies and eigenfunctions, respectively, for vanishing magnetic interaction, $`J_{\mathrm{pd}}=0`$. These have been obtained in the previous section—note that $`H`$ contains the full Coulomb disorder potential. $`\xi _n`$ includes the hole chemical potential.

We rewrite Eq. (4) as

$$H=\underset{n\sigma }{}\xi _nc_{n\sigma }^{}c_{n\sigma }\frac{J_{\mathrm{pd}}}{4}\underset{i}{}\left[(𝐬_i+𝐒_i)^2(𝐬_i𝐒_i)^2\right]$$
(6)
and employ the imaginary-time functional integral. The two spin-squared terms are decoupled using Hubbard-Stratonovich transformations with auxilliary fields $`𝐡_i^\pm `$ coupling to $`𝐬_i\pm 𝐒_i`$. The auxilliary fields are then transformed to new fields $`𝐡_i𝐡_i^++i𝐡_i^{}`$, $`𝐡_i^{}𝐡_i^+i𝐡_i^{}`$. This leads to the Lagrangian

$$L=\underset{n\sigma }{}c_{n\sigma }^{}(_\tau +\xi _n)c_{n\sigma }+iS\underset{i}{}\dot{\varphi }_i\mathrm{cos}\theta _i+\underset{i}{}\frac{𝐡_i𝐡_i^{}}{J_{\mathrm{pd}}}\underset{i}{}(𝐡_i𝐬_i+𝐡_i^{}𝐒_i).$$
(7)
To obtain the mean-field theory we employ a saddle-point approximation. We assume $`𝐡_i`$ and $`𝐡_i^{}`$ to be constant in time but *not* in space, since we wish to retain the disorder effects. The holes and Mn spins are then easily integrated out. For the Mn spins this is best done in the Hamiltonian formalism. This leads to the grand potential

$$\beta \mathrm{\Omega }_0=\mathrm{Tr}\mathrm{ln}\beta G^1(i\omega )\underset{i}{}\mathrm{ln}\frac{\mathrm{sinh}[\beta h_i^{}(S+1/2)]}{\mathrm{sinh}(\beta h_i^{}/2)}+\frac{\beta }{J_{\mathrm{pd}}}\underset{i}{}𝐡_i𝐡_i^{}$$
(8)
with the mean-field hole Green function

$$G_{n\sigma n^{}\sigma ^{}}^1(i\omega )=(i\omega +\xi _n)\delta _{nn^{}}\delta _{\sigma \sigma ^{}}\underset{i}{}\psi _n^{}(𝐑_i)\frac{𝝈_{\sigma \sigma ^{}}}{2}𝐡_i\psi _n^{}(𝐑_i).$$
(9)
The saddle-point equations are obtained by setting the derivatives of $`\beta \mathrm{\Omega }_0`$ with respect to $`𝐡_i`$ and $`𝐡_i^{}`$ to zero,

$`𝐡_i`$ $`=`$ $`J_{\mathrm{pd}}𝐌_i=J_{\mathrm{pd}}{\displaystyle \frac{𝐡_i^{}}{h_i^{}}}SB_S(\beta Sh_i^{}),`$ (10)
$`𝐡_i^{}`$ $`=`$ $`J_{\mathrm{pd}}𝐦_i={\displaystyle \frac{J_{\mathrm{pd}}}{\beta }}{\displaystyle \underset{i\omega }{}}{\displaystyle \underset{n\sigma n^{}\sigma ^{}}{}}G_{n\sigma n^{}\sigma ^{}}(i\omega )\psi _n^{}^{}(𝐑_i){\displaystyle \frac{𝝈_{\sigma ^{}\sigma }}{2}}\psi _n(𝐑_i)`$ (11)
with the usual Brillouin function $`B_S(x)`$. To actually evaluate Eq. (11) we diagonalize the hole sector of the Hamiltonian or, equivalently, the inverse Green function $`G^1`$, writing $`G^1=𝒪g^1𝒪^{}`$, where $`g_{n\sigma }^1=i\omega +\stackrel{~}{\xi }_{n\sigma }`$ is diagonal. We assume the the mean-field Mn spin polarisation $`𝐌_i`$ to be collinear and choose the $`z`$ axis along $`𝐌_i`$. Then $`G^1`$ is symmetric and $`𝒪`$ is orthogonal. Also, $`G^1`$ is then already diagonal in spin space so that we can diagonalize separately in each spin sector. We denote the components of $`𝒪`$ for spin sector $`\sigma `$ by $`𝒪_{nn^{};\sigma }`$. The inversion of $`G`$ is now trivial, $`G=𝒪g𝒪^{}`$. The frequency sum can finally be performed, since $`𝒪`$ does not depend on frequency. Introducing the new eigenfunctions

$$\stackrel{~}{\psi }_{n\sigma }(𝐑_i)=\underset{n^{}}{}\psi _n^{}(𝐑_i)𝒪_{n^{}n;\sigma },$$
(12)
Equation (11) then obtains the form

$$𝐡_i^{}=J_{\mathrm{pd}}\underset{n\sigma }{}n_F(\stackrel{~}{\xi }_{n\sigma })\stackrel{~}{\psi }_{n\sigma }^{}(𝐑_i)\frac{𝝈_{\sigma \sigma }}{2}\stackrel{~}{\psi }_{n\sigma }(𝐑_i).$$
(13)
The mean-field solution is obtained by iterating Eqs. (10) and (13).

Figure 5 shows the magnetisation curves for Mn and hole spins, again for $`x=0.05`$ and $`p=0.3`$ for random and equilibrium configurations. The hole-spin contribution to the total magnetisation is found to be very small. The shape of the Mn magnetisation curve becomes more “normal” and concave from below with annealing , except for a tail at higher temperatures, in good agreement with experimental annealing studies . The anomalous shape for random defects is due to the localisation tendency of the holes, which leads to a shorter-range effective Mn spin interaction. This, together with the random positions of spins, leads to a broad distribution of effective fields acting on these spins. For higher temperatures, only a few spins with strong interactions carry most of the magnetisation, leading to the drawn out tail for random defects. The noise apparent on the curves is due to there being several saddle-points with similar free energy and total magnetisation. We also plot the standard deviations $`\mathrm{\Delta }M`$ of the Mn spin polarisations. The standard deviation becomes comparable with the average Mn spin polarisation at intermediate temperatures below the Curie temperature $`T_c`$. This shows that the reduction of the total magnetisation with increasing temperature is initially dominated by disordering of large local moments. Only at higher temperatures the local moments also decrease.

The Curie temperatures are significantly larger than what one obtains from a virtual crystal approximation (VCA), neglecting disorder completely, see Fig. 5. It is an important open question whether this is an artifact of mean-field theory or indeed a physical result. However, the tails of the magnetisation curves with small average Mn spin polarisation but large standard deviation are dominated by a few anomalously strongly coupled spins. These are overemphasized by mean-field theory (remember that mean-field theory gives an incorrect finite polarisation even for two coupled spins). Thermal fluctuations should be effective in destroying the long-range order in this case, thus reducing $`T_c`$.

Finally, Fig. 6 shows magnetisation curves for equilibrated configurations at various Mn concentrations with the number of holes chosen as above. We see that the curves for smaller Mn and hole concentrations become more similar to the random limit in Fig. 5, which is easily understood since the holes become more localized.

## 5 Conclusions

We have shown that strong Coulomb interactions in highly compensated diluted ferromagnetic semiconductors naturally lead to strong correlations in the positions of charged defects. We also found that these correlations are essential for the description of transport properties, in particular the metal-insulator transition, the persistence of the energy gap, and the shape of the magnetisation curves. It is a pleasure to acknowledge discussions with P. J. Jensen, J. König, U. Krey, W. Nolting, M. E. Raikh, F. Schäfer, and J. Schliemann.

## References