text:
# Untitled Document

Observation of the Charge Symmetry Breaking

d+d4He+π0 Reaction Near Threshold

E.J. Stephenson,1 A.D. Bacher,1,2 C.E. Allgower,1 A. Gårdestig,3

C. Lavelle,1 G.A. Miller,4 H. Nann,1,2 J. Olmsted,1 P.V. Pancella,5 M.A. Pickar,6

J. Rapaport,7 T. Rinckel,1 A. Smith,8 H.M. Spinka,9 and U. van Kolck10,11

1Indiana University Cyclotron Facility, Bloomington, IN 47408

2Department of Physics, Indiana University, Bloomington, IN 47405

3Indiana University Nuclear Theory Center, Bloomington, IN 47408

4Department of Physics, University of Washington, Seattle, WA 98195

5Physics Department, Western Michigan University, Kalamazoo, MI 49008

6Department of Physics and Astronomy, Minnesota State University

at Mankato, Mankato, MN 56001

7Department of Physics and Astronomy, Ohio University, Athens, OH 45701

8Physics Department, Hillsdale College, Hillsdale, MI 49242

9Argonne National Laboratory, Argonne, IL 60439

10Department of Physics, University of Arizona, Tucson, AZ 85721

11RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973

ABSTRACT: We report the first observation of the charge symmetry breaking $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ reaction near threshold at the Indiana University Cyclotron Facility. Kinematic reconstruction permitted the separation of $`{}_{}{}^{4}\mathrm{He}+\pi ^0`$ events from double radiative capture $`{}_{}{}^{4}\mathrm{He}+\gamma +\gamma `$ events. We measured total cross sections for neutral pion production of $`12.7\pm 2.2`$ pb at 228.5 MeV and $`15.1\pm 3.1`$ pb at 231.8 MeV. The uncertainty is dominated by statistical errors.

Charge symmetry is the approximate symmetry of the strong interaction under a specific isospin rotation that interchanges down and up quarks . This symmetry is broken by the different masses of the down and up quarks ($`m_d>m_u`$) and by their electromagnetic interactions. The combination of these two mechanisms leads, for example, to the neutron being heavier than the proton. Within the framework of chiral effective field theory , additional experimental information on the relative contributions of these two mechanisms to charge symmetry breaking (CSB) must, in leading order, come from pion-nucleon scattering. Direct experimental evidence is restricted to elastic scattering and charge exchange experiments with low-energy charged pions where the interpretation is complicated by corrections for the neutron-proton mass difference and electromagnetic interactions between the pions and nucleons . Reactions in which a $`\pi ^0`$ is emitted after being created by one nucleon and rescattered by a second are particularly clean. One example of such a CSB process is the measurement of a forward-backward asymmetry in the cross section for the $`\mathrm{n}+\mathrm{p}\mathrm{d}+\pi ^0`$ reaction .

The $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ reaction violates charge symmetry because the $`\pi ^0`$, whose wavefunction is odd under the interchange of down and up quarks, should not be produced from charge symmetry even (in this case self-conjugate) nuclear systems. The amplitude for CSB is weaker than a similar charge symmetry conserving amplitude by about 1/300, roughly the ratio of the quark mass difference to the nucleon mass. This suggests a $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ total cross section, which depends only on a CSB amplitude, that is as small as tens of picobarns.

Several searches for the $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ reaction have produced only upper limits . A positive report at a deuteron energy of 1.1 GeV  has been questioned because the experiment did not clearly distinguish the photons from $`\pi ^0`$ decay from photons that could have been produced by the double radiative capture $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\gamma +\gamma `$ process. This process is isospin allowed and predicted to appear at the reported cross section . A further search that can clearly distinguish between these two reactions is therefore warranted.

We chose to look for the $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ reaction just above its threshold at 225.5 MeV to avoid other pion producing channels and to take advantage of the clean experimental conditions afforded by the Indiana University Cyclotron Facility’s electron-cooled storage ring. A $`6^{}`$ bend located in one section of the ring provided a site where 4He nuclei, produced in a narrow forward cone just above threshold, could be separated from the circulating deuteron beam. By placing a gas jet target sufficiently upstream of the $`6^{}`$ magnet, it became possible to cover a large solid angle with two arrays of Pb-glass detectors that would be selectively sensitive to photons from the target region.

Figure 1: Layout of the experimental setup showing the target, approximate locations of the Pb-glass arrays, and the magnetic channel in relation to a segment of the electron-cooled storage ring. Quadrupole magnets (Q1, Q2, and Q3), wire chambers (WC1, WC2, and WC3), and scintillation trigger ($`\mathrm{\Delta }`$E1, $`\mathrm{\Delta }`$E2, and E) and veto (V1 and V2) detectors are shown. The luminosity detectors are small and consequently are omitted here.

The layout of the experiment is shown in Fig. 1, with the major features being two arrays of Pb-glass detectors stacked on the left and right sides of the gas target box, the $`6^{}`$ separation magnet, and a magnetic channel consisting of a septum magnet to steer the 4He nuclei away from the downstream ring quadrupole magnets and a quadrupole triplet to confine the 4He nuclei within the acceptance of the detectors at the end of a long, evacuated drift length.

For some $`{}_{}{}^{4}\mathrm{He}+\pi ^0`$ events, the 4He nucleus and both $`\pi ^0`$ decay photons could be recorded. Parts of the solid angle above and below the target were required for differential pumping of the jet target volume. Kinematic constraints just above threshold require very similar photon opening angles and energies for $`{}_{}{}^{4}\mathrm{He}+\pi ^0`$ and $`{}_{}{}^{4}\mathrm{He}+\gamma +\gamma `$ events. Separation of these two reactions relied on measuring the 4He momentum vector using time of flight in the channel for the longitudinal component and scattering angle for the transverse component. The flight time was measured from the first plastic scintillation detector ($`\mathrm{\Delta }`$E1) to the second ($`\mathrm{\Delta }`$E2) at the end of the channel 5.73 m away. A third scintillator (E) in which the 4He nuclei stopped was required for an event trigger along with the absence of any signal in two additional scintillators (V1 and V2). Pulse-height correlations among the trigger scintillators cleanly separated 4He events. Pb-glass information was not used in the trigger. The scattering angle of the 4He was recorded by a multiwire proportional chamber (WC1). Two additional chambers (WC2 and WC3) tracked each 4He nucleus through the channel. From the 4He momentum and the assumption of a two-body final state, it is possible to calculate the missing mass. Double radiative capture events produce a broad distribution of mass values up to a kinematic limit that depends on the beam energy rather than a peak at the $`\pi ^0`$ mass.

The isospin-allowed $`\mathrm{p}+\mathrm{d}{}_{}{}^{3}\mathrm{He}+\pi ^0`$ reaction was used to commission the detector system. The calculation of 3He (and later 4He) momentum from channel time of flight used a model to describe the energy loss (scaled from Janni ) of the 3He passing through the channel detectors and vacuum windows. The time offsets for each photomultiplier associated with the $`\mathrm{\Delta }`$E1 and $`\mathrm{\Delta }`$E2 detectors were adjusted empirically during the subsequent data analysis. For 3He energies close to threshold and runs lasting several hours, a $`\pi ^0`$ mass resolution as small as $`\mathrm{FWHM}=240`$ keV was obtained.

For $`\mathrm{p}+\mathrm{d}{}_{}{}^{3}\mathrm{He}+\pi ^0`$ events in which both $`\pi ^0`$ photons were recorded, photon energies were taken to be the sum of Pb-glass energies starting with the detector recording the highest energy deposition and including its eight nearest neighbors. With summed energies above a threshold chosen to remove most random events and with Pb-glass timing in the correct range relative to the channel trigger, the measured efficiency for detecting the two $`\pi ^0`$ decay photons was $`0.360\pm 0.001`$. Simulations using the GEANT Monte-Carlo program library  and the known $`\pi ^0`$ angular distribution  reproduced the Pb-glass spectral energy shape and agreed with the measured efficiency to $`\pm 0.01`$ This simulation was used to obtain the changes in Pb-glass efficiency between the $`\mathrm{p}+\mathrm{d}{}_{}{}^{3}\mathrm{He}+\pi ^0`$ commissioning run and production running for $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ due to alterations of the kinematics and center-of-mass angular distribution (assuming the $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ cross section to be isotropic). Gains for the Pb-glass detectors were calibrated using cosmic ray muons. This calibration operated continuously during commissioning and production running.

The size of the 3He opening angle cone for $`\mathrm{p}+\mathrm{d}{}_{}{}^{3}\mathrm{He}+\pi ^0`$ calibrated the incident proton energy. Several determinations at slightly different energies were combined by using the frequency of the RF voltage that maintained the beam bunching to calculate the circumference of the storage ring, yielding a value of $`86.786\pm 0.003`$ m that was consistent over time. This calibration was used subsequently for the deuteron beam energy.

The deuterium target was made by directing deuterium gas from a glass nozzle cooled to about 40 K across the circulating storage ring beam ($`2`$ mm wide). The gas flow was adjusted so that the beam lifetime was close to 100 s with the RF voltage and electron cooling on, a value that gave the largest data production rate. The circulating beam current was as large as 2 mA.

The luminosity (product of beam flux and intercepted target thickness) was determined through a separate calibrated monitoring system. For this we used $`\mathrm{d}+\mathrm{d}`$ elastic scattering in the vicinity of $`\theta _{c.m.}=90^{}`$ by placing two scintillator detectors at $`44^{}`$ on either side of the beam. The product of $`\mathrm{d}+\mathrm{d}`$ differential cross section and solid angle for this system was measured in a separate experiment using a molecular HD target. Additional scintillators were added to observe $`\mathrm{d}+\mathrm{p}`$ elastic scattering. Using known $`\mathrm{d}+\mathrm{p}`$ cross sections and the equality of $`\mathrm{d}+\mathrm{p}`$ and $`\mathrm{d}+\mathrm{d}`$ luminosities for an HD target, we obtained the $`\mathrm{d}+\mathrm{d}`$ calibration. Cross sections for $`\mathrm{d}+\mathrm{p}`$ elastic scattering in our energy range have recently become available  with an absolute normalization error of 5.5%. The scintillation detectors provided particle identification information to separate breakup and other backgrounds. An additional position-sensitive silicon detector used to record recoil nuclei from small-angle scattering provided a measurement of the jet target profile as it intercepted the beam. This information was used to make corrections to the scintillator acceptance geometry for both the $`\mathrm{d}+\mathrm{p}`$ and $`\mathrm{d}+\mathrm{d}`$ elastic scattering processes, and to track changes between calibration and production running. The average luminosity during production was $`2.9\times 10^{31}`$ /cm2/s.

Candidate $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ events were required to have the correct pulse height in the three trigger scintillators ($`\mathrm{\Delta }`$E1, $`\mathrm{\Delta }`$E2, and E), usable wire chamber information, and photon signals in both the left and right Pb-glass arrays with an energy sum above threshold and timing coincident with 4He events in the channel. For each candidate event the missing mass was calculated. The results were accumulated into a spectrum, as shown in Fig. 2. Investigations of the quality of the data set showed that random background was essentially removed ($`1`$ event/spectrum); thus we should interpret these results as arising only from $`\pi ^0`$ production and double radiative capture.

Figure 2: Historgrams of the candidate events at the two deuteron bombarding energies as a function of their missing mass value. The smooth curves show the reproduction of these histograms with a Gaussian peak and a continuum.

Production running started at 228.5 MeV, an energy chosen because the half angle of the 4He forward cone ($`1.2^{}`$) was expected to fit safely inside the channel acceptance. Over the much longer $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ production running times, drifts in the time-of-flight measurement became more of a problem, despite efforts to track timing changes by monitoring other particles going along the channel. The width of the missing mass peak was $`\mathrm{FWHM}=510`$ keV, as seen in the top panel of Fig. 2. This meant that the peak was not well separated from the $`{}_{}{}^{4}\mathrm{He}+\gamma +\gamma `$ kinematic upper limit of 136.4 MeV. Below the peak, the flux we would attribute to double radiative capture is attenuated by the acceptance boundaries of the channel, increasing the ambiguity of its shape. The decision was made to complete production running at the higher energy of 231.8 MeV ($`1.75^{}`$ cone) with a kinematic endpoint of 138.0 MeV even though an additional 10% of the $`\pi ^0`$ events would be lost in the channel. The result shown in the lower panel of Fig. 2 still has a broad peak ($`\mathrm{FWHM}=660`$ keV) but there is now a clear distribution of double radiative capture events on either side.

The spectra of Fig. 2 were modelled with a Gaussian peak and a continuum. The continuum shape was obtained from the distribution in missing mass of all 4He events (rate about $`10^3`$ higher than for candidate events alone) reduced by the ratio of the calculated double radiative capture cross section to the phase space value. These 4He events without coincident photons, which may have originated from (d,4He) reactions on residual gas and storage ring structures, were broadly distributed in energy and angle, and were used as an estimate of the product of phase space and channel acceptance. The continuum curves in Fig. 2 are a smooth representation of this shape. The centroids of the Gaussian peaks have a fitting error of less than 60 keV and are consistent with the $`\pi ^0`$ mass. Events in the peak appear to be isotropically distributed in the center of mass as shown by their distribution in scattering angle and time of flight.

The 66 and 50 $`{}_{}{}^{4}\mathrm{He}+\pi ^0`$ events recorded at the two energies of 228.5 and 231.8 MeV lead to total cross section values of $`12.7\pm 2.2`$ pb and $`15.1\pm 3.1`$ pb respectively, including a 6.6% normalization error for all systematic effects. Corrections at the two energies included Pb-glass efficiency (0.34 and 0.32), trigger losses (0.94 and 0.96 for random vetoes), system livetime (0.95 and 0.94), wire chamber efficiency (0.93 and 0.95), and other channel losses from acceptance and multiple scattering (0.95 and 0.81). These cross sections are consistent with being proportional to $`\eta =p_\pi /m_\pi `$ with a combined slope of $`\sigma _{\mathrm{TOT}}/\eta =80\pm 11`$ pb. The integral of the double radiative capture process for the 2 MeV range just below the kinematic limit is $`6.9\pm 0.9`$ pb and $`9.5\pm 1.4`$ pb, values that are more than double predictions .

These results provide the first unambiguous measurement of the $`\mathrm{d}+\mathrm{d}{}_{}{}^{4}\mathrm{He}+\pi ^0`$ cross section as well as that for the competing double radiative capture process. A detailed interpretation of our results will require a careful consideration of different CSB reaction mechanisms. On the basis of chiral power counting arguments, the leading terms are expected to be  those related to $`\pi \eta `$ mixing  and CSB $`\pi +\mathrm{N}`$ scattering . Preliminary plane wave calculations  suggest that $`\pi \eta `$ mixing may be important because all four nucleons contribute coherently, especially if the $`\eta `$ is produced by two nucleons via the exchange of a heavy meson, as in $`\mathrm{p}+\mathrm{p}\mathrm{p}+\mathrm{p}+\pi ^0`$ . The effects of pion scattering are also expected to be enhanced by the initial-state $`\mathrm{d}+\mathrm{d}`$ interactions . Preliminary estimates show that electromagnetic effects are very small.

Both $`\pi \eta `$ mixing and CSB $`\pi +\mathrm{N}`$ scattering are also important in producing an asymmetry about $`90^{}`$ in the $`\mathrm{n}+\mathrm{p}\mathrm{d}+\pi ^0`$ reaction , with the latter determining the sign of the asymmetry. Calculations now underway should allow the extraction of both processes from the combined experimental results. Since both depend on the quark mass difference, experimental information that separates each contribution will greatly enhance efforts to relate the hadronic matrix elements to the underlying QCD theory.

We would like to thank V. Anferov, G.P.A. Berg, B. Chujko, J. Doskow, G. East, C.C. Foster, W. Fox, D. Friesel, T. Hall, C.J. Horowitz, A. Kuznetsov, V. Medvedev, H.O. Meyer, D. Patalahka, R.E. Pollock, A. Prudkoglyad, S. Shastry, T. Sloan, K. Solberg, J. Vanderwerp, and B. von Przewoski for their contributions to the preparation of this experiment. Our work was supported in part by the National Science Foundation grants PHY-97-22538, PHY-98-02872, and PHY-00-70368, the Department of Energy grant DE-AC02-98CH10886, the Department of Energy OJI award DE-FG03-01ER41196, RIKEN, Brookhaven National Laboratory, and an Alfred P. Sloan Fellowship.

1. G.A. Miller, B.M.K. Nefkens, and I. Šlaus, Phys. Rep. 194, 1 (1990).

2. G.A. Miller and W.T.H. van Oers, in Symmetries and Fundamental Interactions in Nuclei, eds. W.C. Haxton and E.M. Henley (World Scientific, Singapore, 1995) p. 127; and references therein.

3. U. van Kolck, Few Body Systems, Suppl. 9, 444 (1995).

4. U. van Kolck, J.A. Niskanen, and G.A. Miller, Phys. Lett. B 493, 65 (2000).

5. W.R. Gibbs, Li Ai, and W.B. Kaufmann, Phys. Rev. Lett. 74, 3740 (1995).

6. C. Matsinos, Phys. Rev. C 56, 3014 (1997).

7. A.K. Opper and E. Korkmaz, TRIUMF experiment 704 proposal.

8. J. Banaigs et al., Phys. Rev. Lett. 58, 1922 (1987); and references therein.

9. L. Goldzahl, J. Banaigs, J. Berger, F.L. Fabbri, J. Hüfner, and L. Satta, Nucl. Phys. A533, 675 (1991).

10. D. Dobrokhotov, G. Fäldt, A. Gårdestig, and C. Wilkin, Phys. Rev. Lett. 83, 5246 (1999).

11. J.F. Janni, At. Data and Nucl. Data Tables 27, 150 – 529 (1982).

12. GEANT program library, http://wwwinfo.cern.ch/asdoc/geant.html3/geantall.html.

13. M.A. Pickar et al., Phys. Rev. C 46, 397 (1992).

14. K. Ermisch, Ph.D. thesis, Rijkuniversiteit Groningen, 2003.

15. A. Gårdestig, private communication.

16. A. Fonseca, A. Gårdestig, C. Hanhart, C.J. Horowitz, G.A. Miller, A. Nogga, and U. van Kolck, private communication.

17. S.A. Coon and B.M. Preedom, Phys. Rev. C 33, 605 (1986).

18. T.-S.H. Lee and D.O. Riska, Phys. Rev. Lett. 70, 2237 (1993).


format: text

text:
# Power-law versus exponential distributions of animal group sizes

## 1 Introduction

It is hard to find animals in nature that do not aggregate (school, herd, swarm, or flock). Despite the universality of aggregation and its ecological and economic importance for the estimation of wildlife abundance, the statistical question about animal group-size has been involved in much confusion. Theoretically the Gibbs-Boltzmann (exponential) distribution has been proposed for the animal group-size distribution by applying a maximum entropy principle (Okubo, 1986), whereas practically the lognormal distribution has been used in fisheries (Matsuishi et al., 1993). Only lately, it has been suggested that power-law distributions may be quite generic (Bonabeau & Dagorn, 1995). Anderson (1981), a pioneer in statistical research into group formation, proposed a stochastic dynamic equation for the size of fish school, in which the possibility of power laws was already presented but not exploited. Bonabeau & Dagorn (1995) presented a model of animal aggregation inspired by a physical model of particle aggregation (Takayasu, 1989), and predicted that the group-size distributions follow a power-law decay with exponent 3/2. The empirically determined power-law size distributions of fish schools have exponents in the range from 0.7 to 1.8 (Niwa, 1998; Bonabeau et al., 1999).

$`N`$ animals moving together form an $`N`$-sized group. The group-size distribution $`W(N)`$ is proportional to the observed number of $`N`$-sized groups. Power-law distributions, $`W(N)N^\beta `$, have somewhat unusual properties. They do not have a well-defined mean when $`\beta 2`$. Fat-tailed group-size distributions are necessarily truncated at a cut-off size because the population is finite, but truncated power law must be distinguished from purely rapidly decreasing ones, as they exhibit specific properties, e.g., a violation of the central-limit theorem (Mantegna & Stanley, 1995; Axtell, 2001). In order to investigate the possibility of power-law regimen in the group-size distribution, let us reconsider some existing data (Table 1) in terms of the population distribution $`P(N)=NW(N)`$, which is proportional to the ratio of animal population in $`N`$-sized groups to total population.

Suppose that we have a data set arranged in a frequency distribution with $`n`$ classes. Let $`W_i\mathrm{\Delta }N`$ be the frequency of animal groups observed in the $`i`$-th class (where $`\mathrm{\Delta }N`$ is the class width); the class mark of group size, $`N_i`$, is the midpoint of the class. $`W_i`$ then represents the distribution density of group sizes. The number of individuals associated with groups of the $`i`$-th class is given by $`N_iW_i\mathrm{\Delta }N`$ (denoted by $`P_i\mathrm{\Delta }N`$). $`P_i`$ then represents the distribution density of population. Shown in Fig.1 is a semilogarithmic plot of the population distribution of pelagic fishes. Surprisingly, I find that the scaled population distributions appear to collapse onto a single curve, an exponentially decreasing function. It implies that (i) the population distribution has a well-defined mean

$$N_P=\underset{i=1}{\overset{n}{}}N_iP_i/\underset{i=1}{\overset{n}{}}P_i$$
(1)
which is hereafter referred to as $`P`$-mean; (ii) the school-size distribution displays robust scale-invariant behavior with the power-law index $`\beta =1`$; and (iii) the power-law distribution of school sizes is truncated at the cut-off size equal to $`P`$-mean. $`P`$-mean is the average of the population distribution among group sizes, i.e., the expected size of the groups in which an arbitrary individual engages. The population distribution is the crossover function from power-law to exponential decay.

## 2 The Model

The empirical finding for fish school size suggests that the population distributions among school sizes are identical for fishes regardless of large diversities in behavioral and ecological conditions. I propose a simple stochastic-differential-equation model accounting for group-size statistics. Assume binary splitting independent of group size (breakup rate $`p`$) and coalescing on contingent encounter. Trace the size change of the group that a certain individual (named “A”) rides. The group size at time $`t`$ is denoted by $`N_\text{A}(t)`$. A discrete countable number of group size is replaced by a continuous group size, providing that there is a sufficiently large reservoir of individuals. Rewriting the difference (the size change $`\mathrm{\Delta }N_\text{A}`$ for a finite time interval $`\mathrm{\Delta }t`$) to the differential of the group size for an infinitesimal time interval $`\text{d}t`$,

$$\mathrm{\Delta }N_\text{A}(t)=N_\text{A}(t+\mathrm{\Delta }t)N_\text{A}(t)\text{d}N_\text{A}(t)$$
(2)
gives the Ito stochastic differential equation for group size:

$$\text{d}N_\text{A}(t)=\frac{p}{2}\left(N_\text{A}(t)N_P\right)\text{d}t+\sigma \left(N_\text{A}(t)\right)\text{d}B(t)$$
(3)
where $`\text{d}B`$ is Wiener noise term of normal distribution with zero mean and variance $`\text{d}t`$, and $`\sigma \left(N_\text{A}\right)`$ is the standard deviation of group-size changes. Fluctuations of the group size in aggregation-breakup processes, $`\sigma \left(N_\text{A}\right)`$, are numerically determined as follows:

Assume for simulations that there are $`s`$ sites, coarse-grained zones of space, on which $`\mathrm{\Phi }`$ individuals move. The number of individuals, $`\mathrm{\Phi }`$, is conserved. Each site is either empty or occupied by one group. At each discrete time step, all groups move towards a randomly selected site. They may move to any site with equal probability. This corresponds to the mean-field theory. When $`M`$\- and $`N`$-sized groups happen to move to the same site, they aggregate to form an $`(M+N)`$-sized group. Each group with a size greater than or equal to 2 splits into a pair of groups with the probability $`\widehat{p}`$ at each time step. It is assumed that the probability $`\widehat{p}`$ for a group to split is independent of its size and that the sizes of splitting groups are uniformly distributed: a probability for an $`N`$-sized group to split into $`M`$\- and $`(NM)`$-sized groups is represented by

$$K_\text{b}(N|M,NM)=K_\text{b}(N)=\frac{\widehat{p}}{N1}$$
(4)
for $`N2`$. Then, the expected decrement and increment of size of the group that a certain individual “A” rides are $`\widehat{p}N_\text{A}/2`$ and $`\mathrm{\Phi }/s`$ (i.e., coarse-grained spatial population density) at each time step, respectively. Monte-Carlo simulations of animal-group aggregation show that fluctuations of the group size in aggregation-breakup processes exponentially depend on the size (Fig.2), and give

$$\sigma ^2(N)=2D\mathrm{exp}\left(\frac{N}{N_P}\right)$$
(5)
where $`D`$ denotes the size of fluctuations. The same results can be obtained with another model for a group to split: the size distribution of splitting groups is binomial.

$$K_\text{b}(N|M,NM)=\frac{\widehat{p}}{2^N2}\left(\genfrac{}{}{0pt}{}{N}{M}\right)$$
(6)
Let the domain of the size $`N`$ be real numbers. The aggregation-breakup dynamics \[eqn (3)\] is supposed to be symmetric about $`N=0`$. By performing the change of variable (Richmond, 2001)

$$\frac{\text{d}x}{\text{d}N}=\frac{\sqrt{p/2}}{\sigma (N)}$$
(7)
eqn (3) becomes the Langevin equation for a Brownian particle moving in a potential $`U(x)`$:

$$\text{d}x(t)=\frac{U}{x}\text{d}t+\sqrt{\frac{p}{2}}\text{d}B(t)$$
(8)
where

$$U(x)=\frac{p}{4}\mathrm{ln}\left(1\frac{x}{N_P\sqrt{p/D}}\right)\left[1\frac{pN_P^2}{D}\left(1\frac{x}{N_P\sqrt{p/D}}\right)^2\right]$$
(9)
The induced potential is schematically depicted in Fig.3. If the potential $`U(x)`$ is continuously differentiable, the following relation is obtained:

$$pN_P^2=D$$
(10)
Therefore the breakup rate $`p`$, which corresponds to the dissipation in the stochastic differential equation, and the noise $`D`$ cannot be independent. The numerical results indicate that this fluctuation-dissipation relation is already satisfied (Fig.2).

## 3 Results

The stationary solution of eqn (3) reads the probability for individual “A” to be found in $`N`$-sized groups, which is proportional to the stationary population distribution among group sizes, $`P(N)`$. The distribution $`P(N)`$ has the desired form and the stationary group-size distribution $`W(N)[=N^1P(N)]`$ is rigorously written as

$$W(N)N^1\mathrm{exp}\left[\frac{N}{N_P}\left(1\frac{e^{N/N_P}}{2}\right)\right]$$
(11)
where eqn (5) is adopted for the variance with eqn (10). There is no parameter to be adjusted, since $`P`$-mean is obtained directly from the data \[because of the fluctuation-dissipation relation, eqn (10), parameters $`p`$ and $`D`$ do not explicitly appear in the group-size distribution\]. The self-consistency in the average of population distribution

$$\frac{{\displaystyle _0^{\mathrm{}}}q\mathrm{exp}\left[q\left(1e^q/2\right)\right]\text{d}q}{{\displaystyle _0^{\mathrm{}}}\mathrm{exp}\left[q\left(1e^q/2\right)\right]\text{d}q}=1$$
(12)
is numerically certified by using $`2^8`$-digit precision for internal computations, where $`q=N/N_P`$.

The model accounting for group-size statistics presented here is consistent with a number of numerical experiments. Fig.4 represents the group-size distribution obtained from simulations. The model is tested with empirical data for the fish school-size distributions and a close agreement is found (Fig.5). The normalization factor is given by

$$\left[_0^{\mathrm{}}\mathrm{exp}\left[q\left(1e^q/2\right)\right]\text{d}q\right]^1=0.881237$$
(13)
The size distribution $`W(N)`$ follows a power law with exponent $`\beta =1`$ to a cut-off size $`N_P`$, i.e., a well-defined mean of the associated distribution $`P(N)`$ which is exponentially decaying function. Naturally, a well-defined mean does not exist for the group-size distribution. Moreover, a coarse-grained description of the group aggregation suggests that $`P`$-mean varies in proportion as the number of individuals or the spatial population density $`\rho `$, and inversely as the breakup rate:

$$N_P\frac{\rho }{p}$$
(14)
which is numerically approved (Fig.6). The group-size distribution, therefore, follows a universal scaling law, i.e., it decays as a power law with an exponential truncation controlled not only by the system size (number of individuals) but also the breakup rate. I indeed find a collapse of the empirical data onto a single curve as well as the collapse of the numerical results. Eqn (14) is utilized for data normalization (Table 1).

## 4 Discussion

I have reconsidered some existing data on pelagic fishes and found identical statistical signatures: the school-size distribution follows a power-law decay with the exponent 1 up to a cut-off size which is given by $`P`$-mean. The rationale behind this robust scaling invariance is the exponential fluctuation of school-size change. I have proposed an Ito stochastic differential equation governing the evolution of the group size that a certain individual rides, which is quite different approach from Anderson’s (1981) model based on a stochastic dynamics for the change of the number of individuals in a certain group. The model predicts not only the power-law behavior observed in nature, but also the deviation from pure power-law towards exponential decay. The similarity between the empirical distribution and the distribution obtained from the model is striking, though there is no fitting parameter. The exponential fluctuation of group size is the essential ingredient of underlying aggregate-breakup dynamics that influence animal group size. If the size fluctuation $`\sigma ^2(N)`$ were proportional to the group size $`N`$, the group-size distribution would follow a power-law decay with exponent 3/2, as was suggested by Bonabeau & Dagorn (1995).

The interacting group systems introduced here have a marked feature in contrast to physical systems in thermal equilibrium: the aggregate-breakup process in animal-group systems does not allow for detailed balance as below; while Okubo (1986) made the detailed balance assumptions in order to obtain the group-size distribution.

The Smoluchowski rate equation is the alternative equation governing the group-size distribution (Gueron & Levin, 1995; Gueron, 1998; Durrett et al., 1999), and is equivalent to the following master equation for associated distribution of population (Niwa, 1998):

$`P(N,t+\mathrm{\Delta }t)P(N,t)`$ (15)
$`=P(N,t){\displaystyle \underset{M=1}{\overset{\mathrm{\Phi }}{}}}m(N,M;t)+{\displaystyle \underset{M=1}{\overset{\mathrm{\Phi }}{}}}P(M,t)m(M,N;t)`$
with the transition probability for a specific individual which engages in a group of size $`N`$ to pass into a group of size $`M`$ in a finite time interval $`\mathrm{\Delta }t`$:

$`m(N,M;t)`$ (16)
$`=K_\text{b}(N|M,NM){\displaystyle \frac{M}{N}}+P(MN,t){\displaystyle \frac{K_\text{a}(N,MN|M)}{MN}}`$
where $`K_\text{a}`$ denotes the probability for an $`N`$-sized group to join with an $`M`$-sized group in a finite time interval $`\mathrm{\Delta }t`$. Whenever two groups meet they are supposed to join, which implies $`K_\text{a}=s^1`$ for a coarse-grained model of animal-group system. Assume the detailed balance condition: there are as many transitions per $`\mathrm{\Delta }t`$ from $`M`$-sized groups to $`N`$-sized groups as from $`N`$ to $`M`$ by the inverse process. The following recurrence formula for stationary solution is then obtained \[$`K_\text{b}`$ is given by eqn (4)\]

$$P(N)=P(N1)N\frac{P(1)}{s\widehat{p}}$$
(17)
which leads to an explicit solution for group-size distribution

$$W(N)=\frac{W(1)}{e^\mu \pi _{N1,\mu }}$$
(18)
where $`\pi _{N,\mu }`$ is the Poisson distribution with mean $`\mu =s\widehat{p}/W(1)`$. The result contradicts the simulated group-size distribution; besides, eqn (18) violates the conservation law of population. The detailed balance conditions therefore do not hold in the aggregate-breakup systems, or else we may have another deterministic equation which is consistent with micro-reversibility.

The model developed here can apply to a wide spectrum of cases on animal species. Fig.7 shows the herd-size distributions for mammalian herbivores; these distributions are much better fitted by a power law than by an exponential decay, as was suggested by Okubo (1986). We may however notice that the empirical data are not in perfect agreement with eqn (11). It may be explained by biases in the data. Another possibility bringing a subtle difference between data and theory is that the uniform-breakup assumption (splitting probability independent of group size) is no longer valid.

There is yet another possibility. This subtle difference may be caused by the fluctuation-dissipation relation reduced above. It can be generalized as follows:

$$pN_P^2=x_0^2D$$
(19)
where $`x_0^21`$, and potential $`U(x)`$ is no longer differentiable at $`x=0`$ but continuous (Fig.3). The population distribution among group sizes is then modified into

$$P(N)\mathrm{exp}\left[\frac{N}{N_P}\left(1\frac{pN_P^2}{2D}e^{N/N_P}\right)\right]$$
(20)
which guarantees the self-consistency in the average of distribution as well as eqn (12). The modified size distribution by adopting eqn (19) matches better the data. Yet, the modified distribution \[eqn (20)\] can apply to schools of sardinellas caught in the up-welling areas (Bonabeau et al., 1999), and the agreement between the empirical distribution and the model’s prediction is remarkable (Fig.8).

Figure legends

* Empirical population distribution of pelagic fishes, as shown in Table 1 (except sardinellas). The scaled distributions $`P_iN_P`$ are plotted against the scaled school sizes $`N_i/N_P`$. The data are clearly fitted by an exponential decay, suggesting that the distributions are identical, e.g., $`\mathrm{exp}\left(N/N_P\right)`$ (broken line). The proposed model is somewhat more complicated (solid line).
* Exponential fluctuations of group-size change (numerical results). The ordinate is the variance of size change of the group, $`\left(\mathrm{\Delta }N_\text{A}\mathrm{\Delta }N_\text{A}_N\right)^2_N`$, for a finite time interval (single time step of simulation run), divided by $`\widehat{p}N_P^2`$, where $`\mathrm{\Delta }N_\text{A}_N`$ denotes the expected value of the size change $`\mathrm{\Delta }N_\text{A}`$ for a finite time interval, providing that individual “A” rides an $`N`$-sized group. The abscissa is the group size $`N`$ divided by $`P`$-mean. Simulations show that the group sizes exponentially fluctuate \[solid line represents $`\sigma ^2(N)\mathrm{exp}\left(N/N_P\right)`$\]. For the simulation parameters, consult Table 2. Note that the variance is scaled based on the fluctuation-dissipation relation \[described later; eqn (10)\].
* Potential induced by the change of variable. (a) Potential $`U(x)`$ is continuously differentiable (solid line). The ordinate is the potential $`U(x)`$ divided by the breakup rate $`p`$. The abscissa is the reduced size $`x`$. $`P`$-mean $`N_P`$ takes positive values for $`x0`$ and negative values for $`x<0`$. The reduced domain of $`x`$ is $`(1,1)`$. If potential $`U(x)`$ is continuous (broken line showing $`U/p`$ versus $`x/\sqrt{pN_P^2/D}`$), the fluctuation-dissipation relation is generalized into eqn (19). (b) The group size $`N`$ as a function of the reduced size $`x`$. The ordinate is the scaled group size $`N/N_P`$.
* Simulated group-size distribution. The scaled distributions $`W_NN_P`$ are plotted against the scaled group sizes $`N/N_P`$. The scaled distributions collapse onto a single curve that corresponds to eqn (11) with normalization factor eqn (13). The simulation parameters are summarized in Table 2.
* Empirical school-size distribution of pelagic fishes (the same data sets as Fig.1). The scaled distributions $`W_iN_P`$ are plotted against the scaled school sizes $`N_i/N_P`$. The scaled data collapse onto a single curve that corresponds to eqn (11) with normalization factor eqn (13).
* Plot of $`P`$-mean against the scaled population size $`\mathrm{\Phi }/s\widehat{p}`$, i.e., the ratio of population density to breakup probability (numerical results). Parameters are summarized in Table 2. The broken line is a fit with slope $`2.907\pm 0.017`$.
* Empirical herd-size distribution of mammalian herbivores, as shown in Table 1 (scaled size distributions). The solid and broken lines correspond to eqn (11) with normalization factor \[eqn (13)\] and to a least-squares fit of the modified size distribution with $`x_0^2=7.96\pm 1.61`$, respectively.
* Two kinds of empirical size distributions of sardinellas schools, as shown in Table 1 (scaled size distributions). (a) Population distribution. A least-squares fit of the data to eqn (20) with $`x_0^2=14.14\pm 1.95`$ is shown. (b) School-size distribution. The solid and broken lines correspond to eqn (11) with normalization factor \[eqn (13)\] and to a fit of the modified school-size distribution, respectively.

Figures

* Empirical population distribution of pelagic fishes
* Exponential fluctuations of group-size change
* Potential induced by the change of variable
* Simulated group-size distribution
* Empirical school-size distribution of pelagic fishes
* Plot of $`P`$-mean against the scaled population size
* Empirical herd-size distribution of mammalian herbivores
* Two kinds of empirical size distributions of sardinellas schools


text:
# A Model-Independent Determination of the Expansion and Acceleration Rates of the Universe as a Function of Redshift and Constraints on Dark Energy

## 1 Introduction

A traditional task of cosmology is to determine the global geometry and dynamics of the universe. The field has been revolutionized by the modern measurements of CMBR fluctuations (e.g., Bennett et al. 2003, Spergel et al. 2003, and references therein), the use of distant supernovæ (SNe) in a Hubble diagram (see, e.g., Riess 2000, Leibundgut 2001, and references therein), radio galaxies (e.g. Daly & Guerra 2002), and many other advances. What these modern measurements have now convincingly demonstrated is that the global mass/energy budget of the universe, and thus its dynamics, is dominated by a so-called dark energy component, which accounts for $`\text{ }>\text{ }70`$% of the closure density today. Einstein’s cosmological constant, $`\mathrm{\Lambda }`$, is one special (and viable) case. More generally, this mysterious dark energy component is characterized through the equation of state, $`w=p/\rho `$, where $`p`$ is the pressure and $`\rho `$ the energy density; the cosmological constant solution corresponds to $`w=1`$. For reviews and further references, see, for example, Sahni & Starobinsky (2000), Turner (2002a, 2002b), Peebles & Ratra (2003), and Padmanabhan (2003).

The nature of the dark energy (including its evolution in redshift, if any) is one of the most outstanding problems of physics and astronomy today. Constraining it through analysis of cosmological data is a task of a critical importance, and every new data set or analysis method can provide valuable insights into this problem.

Several recent studies have focused on the use of supernovae to determine the properties of the dark energy (Starobinsky 1998; Huterer & Turner 1999; Saini et al. 2000; Chiba & Nakamura 2000; Maor, Brustein, & Steinhardt 2001; Goliath et al. 2001; Astier 2001; Gerke & Efstathiou 2002; Weller & Albrecht 2002; and Padmanabhan & Choudhury 2002). The key ingredients are luminosity distances to sources over a broad range of redshift, preferably including sources at high redshift. Most of these analyses have focused on constraints on an evolving scalar field such as that used to define quintessence (Caldwell, Dave, & Steinhardt) or a rolling scalar field (Peebles & Ratra 1988). More recently, other types of models have been proposed to account for the acceleration of the universe, such as stringy dark energy (Frampton 2002) and k-essence (Armendariz-Picon, Damour, & Mukhanov 1999; Barger & Marfatia 2001).

Here, we focus on direct empirical determinations of the dimensionless expansion rate $`E(z)`$ and acceleration rate $`q(z)`$ as functions of redshift. These require values for dimensionless coordinate distances to sources over a broad range of redshifts. We provide both a new data set, and a new method for estimating of $`E(z)`$ and $`q(z)`$.

We first present coordinate distances to 20 radio galaxies (RGs), reaching out to $`z1.8`$, and thus supplementing the existing SN data in what is a critical redshift regime. These RG data can be used to compare model predictions of any flavor of dark energy with the observations. These coordinate distances are derived and listed in §2. For completeness, and to compare the RG and supernova SN results, the coordinate distances to 78 SNe are also listed in §2. In §3, we derive the expressions for a direct determination of $`E(z)`$ and $`q(z)`$ from measurements of the dimensionless coordinate distances $`y(z)`$. In §4, we describe our simple numerical differentiation technique which can be used to implement these concepts on the real data. We illustrate the method and present our preliminary results based on the current RG and SN data sets in §5, and discuss implications for the properties of the dark energy in §6. A summary and discussion follows in §7.

## 2 Dimensionless Coordinate Distances

The values of coordinate distances to sources at high redshift can be used to determine or constrain global cosmological parameters, and to understand the properties and redshift evolution of the dark energy. Coordinate distances $`(a_or)`$ may be obtained from luminosity distances $`d_L`$ or angular size distances $`d_A`$, since these are simply related to the coordinate distance: $`d_L=(a_or)(1+z)`$, and $`d_A=(a_or)/(1+z)`$ (e.g. Weinberg 1972).

The dimensionless coordinate distance, $`y(z)`$ is simply related to the coordinate distance $`a_or`$, $`y(z)=H_0(a_or)`$ (e.g., Carroll, Press, & Turner 1992; Peebles 1993). The luminosity distance $`d_L`$ and the angular size distance $`d_A`$ are also simply related to the dimensionless coordinate distance: $`d_L=H_0^1y(z)(1+z)=H_0^1D_L`$, where $`D_L`$ is the dimensionless luminosity distance (e.g. Perlmutter et al. 1999), and $`d_A=H_0^1y(z)/(1+z)`$.

Observations of type Ia supernovae and type IIb radio galaxies allow estimates of the dimensionless coordinate distances to sources at different redshift.

The use of FRIIb radio galaxies to determine the angular size distance or coordinate distance to radio galaxies at different redshifts is described in detail elsewhere (e.g. Podariu et al. 2003; Daly & Guerra 2002; Guerra, Daly, & Wan 2000; Daly 1994). In addition to the use of FRIIb radio galaxies addressed here, other methods of using radio galaxies and quasars to determine coordinate distances are discussed by Buchalter et al. (1998), Gurvits, Kellermann, & Frey (1999), Vishwakarma (2001), Lima & Alcaniz (2002), and Chen & Ratra (2003). Here FRIIb radio galaxies are used to obtain dimensionless coordinate distances to 20 radio galaxies following the method described, for example, by Daly & Guerra (2002).

In the radio galaxy method proposed by Daly (1994), one model parameter $`\beta `$ enters into the ratio $`R_{}<D>/D_{}`$; this ratio also depends on observed quantities and the dimensionless coordinate distance $`y(z)`$. In this model the ratio $`R_{}`$ is equal to a constant, $`\kappa `$:

$$R_{}(\beta ,y(z))=\kappa .$$
(1)
The constants $`\kappa `$ and $`\beta `$ and their uncertainties are obtained by fitting all of the data to equation (1), as described in detail by Guerra, Daly, & Wan (2000), and Daly & Guerra (2002). The ratio is given by

$$R_{}=k_oy^{(6\beta 1)/7}(k_1y^{4/7}+k_2)^{(\beta /31)},$$
(2)
where $`k_o`$, $`k_1`$, and $`k_2`$ are observed quantities (described in detail in the Appendix of Guerra, Daly, & Wan 2000). Equation (1) with $`R_{}`$ given by equation (2) allows a determination of $`y(z)`$ to each source; $`y(z)`$ is implicitly known for each source and is determined using an iterative technique. The values of $`y`$ obtained along with the one $`\sigma `$ error of y are listed in Table 1. In determining the one $`\sigma `$ error bar on $`y(z)`$, the uncertainties of $`\kappa `$, $`k_o`$, $`k_1`$, and $`\beta `$ have been included; $`k_2`$ is known to high precision as it is the energy density of the microwave background radiation at the source redshift, and is the term that describes the effects of inverse Compton cooling of relativistic electrons by the microwave background radiation. The best fit values of $`\kappa `$ and $`\beta `$ vary slightly depending upon whether just the radio galaxy data are fit, or whether both the radio galaxy and supernova data are fit. Values of $`y`$ obtained using the best fit parameters to radio galaxies alone are labeled $`y_s`$ in Table 1, and those obtained using the best fit to both the radio galaxy and supernova data are labeled $`y_j`$. That is, the best fit values of $`\kappa `$ and $`\beta `$ change slightly depending upon whether just the radio galaxy data are fit (referred to with a subscipt “s”), or whether the radio galaxy and supernovae data are fit simultaneously (referred to with a subscript “j”). Best fit values of $`\kappa `$ and $`\beta `$ are listed in Table 2. Note that the radio galaxy method does not rely upon a low-redshift normalization; the best fit values of $`\kappa `$ and $`\beta `$ are determined using all of the data.

The best fit values of $`\kappa `$, $`\beta `$, and $`_{}`$ (described below) and their error bars are included in Table 2, where the 54 supernovae included in the “primary fit C” of Perlmutter et al. (1999) and the 20 radio galaxies discussed here were studied. Values obtained from the fits of Daly & Guerra (2002) that allow for quintessence in a spatially flat universe with separate (s) and joint (j) fits to the radio galaxy and supernovae data are labelled “Q.” Best fit values obtained in the rolling scalar field model analyzed by Podariu et al. (2003) are labelled “SF.” As the number of data points in the fit increases, the value of each constant, $`\kappa `$, $`\beta `$, and $`_{}`$, becomes independent of the assumptions of the fit. For example, the 54 supernovae points yield a consistent value of $`_{}`$ for fits that include supernovae only, or radio galaxies and supernovae, and in a universe with quintessence or a rolling scalar field. Since the value of $`_{}`$ changes so little when fit in different models, and when fit including or excluding radio galaxies, only one value of $`y`$ is listed for each supernova. The 20 radio galaxy points show some small variations in the values of $`\kappa `$ and $`\beta `$ obtained with radio galaxies alone, or radio galaxies and supernovae, obtained in a universe with quintessence or a rolling scalar field. As more radio galaxy data points are added, the values of the constants will be more accurately determined. New runs were done that include the full 78 supernovae listed here and the best fit parameters and their error bars are the same as those listed in Table 2.

The coordinate distances to the supernovae are determined following the procedures of Perlmutter et al. (1999) and Riess et al. (1998). In the application of supernovae type Ia as a distance indicator, there is one model parameter $`\alpha `$ which is used to determine the effective apparent B band magnitude at maximum brightness $`m_B^{eff}`$. This is related to the dimensionless coordinate distance $`y(z)`$:

$$m_B^{eff}(\alpha )=_{}+5log[(1+z)y(z)].$$
(3)
The constant $`_{}`$ is determined by fitting all of the supernova data, and is simply related to the absolute magnitude of the supernova peak brightness $`M_B`$: $`_{}=M_B+255logH_0`$ (see Perlmutter et al. 1999). Equation (3) is then used to determine $`y(z)`$ to each of the 54 supernovae in the “primary fit C” of Perlmutter et al. (1999), the 37 supernovae presented by Riess et al. (1998), and the 1 high-redshift supernova published by Reiss et al. (2001), with the magnitude of this source corrected for gravitational lensing (Benitez et al. 2002). The one $`\sigma `$ uncertainty of y is obtained by combining the uncertainties of $`_{}`$ and $`m_B^{eff}`$. These values of listed in Table 2. There are 14 sources that are present in both the Riess et al. (1998) and Perlmutter et al. (1999) samples used here. In the determinations of $`E(z)`$ and $`q(z)`$, average values of $`y`$ with appropriate error bars were used for these duplicate sources; these values are listed in Table 4. The values of $`y(z)`$ are shown in Figures 1 and 2. The good agreement between coordinate distances determined using radio galaxies and supernovae at similar redshifts is easy to see in these figures.

To test the reliability of the values of $`y(z)`$ obtained for the radio galaxies, a comparison was made between cosmological parameters obtained directly from the radio galaxy data alone and those obtained from the values of $`y_s`$ listed in Table 1 in a quintessence model (see line 1 of Table 2). Each value of $`y_s`$ was substituted into equation (3) to obtain an equivalent effective apparent magnitude for the radio galaxy; the value of $`_{}`$ obtained for supernovae alone in a universe with quintessence (line 1 of Table 2) was adopted. These effective apparent magnitudes were then analyzed in a universe with quintessence, and the best fit parameters and their one sigma ranges compared with those obtained directly from the radio galaxy data. First, the $`\chi ^2`$ per degree of freedom went from 16.5/16 to 15.6/15, so the reduced $`\chi ^2`$ remains fairly constant; the number of degrees of freedom drops by one in the new fit since one new parameter, $`_{}`$ is fit. The one sigma range of $`\mathrm{\Omega }_m`$ is 0.0 to 0.24 in the original fit, and is 0.0 to 0.17 in the new fit. The one sigma range of $`\mathrm{\Omega }_Q`$ is 0.76 to 1.0 in the original fit, and is 0.83 to 1.0 in the new fit. The one sigma range of $`w`$ in the original fit is -1.3 to -0.43 centered on -0.73, and is -1.5 to -0.56 centered on -0.8 in the new fit. The best fit value of $`_{}`$ is 23.83 $`\pm 0.08`$, compared with the input value $`_{}`$ of 23.91 $`\pm 0.03`$ used to define an effective apparent magnitude for each radio galaxy. Thus, the cosmological parameters obtained directly from the radio galaxies are very similar to those obtained from the values of $`y(z)`$ listed in Table 1.

## 3 Computation of E(z) and q(z) from the Coordinate Distances

The determinations of the dimensionless coordinate distances do not require any assumptions regarding cosmological parameters, the dark energy, or the redshift evolution of these components once the values of the constants $`\kappa `$, $`\beta `$, and $`_{}`$ have been determined. The first and second derivatives of the dimensionless coordinate distance with respect to redshift can be used to construct a model-independent determination of the dimensionless expansion rate $`E(z)=H(z)/H_0`$, and the acceleration rate $`q(z)=\ddot{a}a/\dot{a}^2`$.

These follow from the relation between redshift $`z`$ and the cosmic scale factor $`a(t)`$, $`(a(t)/a_o)=(1+z)^1`$, and the Robertson-Walker line element, which describes a homogeneous isotropic expanding universe,

$$d\tau ^2=dt^2a^2(t)\left(\frac{dr^2}{1kr^2}+r^2d\theta ^2+r^2sin^2\theta d\varphi ^2\right)$$
(4)
(see, for example, Weinberg 1972). It is well known that these imply $`H(z)\dot{a}/a=\sqrt{1kr^2}(dy/dz)^1H_0,`$ and with $`k=0`$ and $`H(z)=H_0E(z)`$,

$$E(z)=(dy/dz)^1$$
(5)
(e.g. Weinberg 1972; Peebles 1993). Recent CMB measurements indicate that our universe has zero space curvature, $`k=0`$ (e.g., Bennett et al. 2003, Spergel et al. 2003). Thus, in principle, the data $`y(z)`$ can be used to empirically determine the the $`dy/dz`$ and the dimensionless expansion rate $`E(z)`$. This, in turn, is related to cosmological parameters, such as dark energy, and their redshift evolution as discussed in §6. For example, in a universe with quintessence (Caldwell, Dave, & Steinhardt 1998), which has a time-independent equation of state $`w=P/\rho `$, $`E^2(z)=\mathrm{\Omega }_i(1+z)^{n_i}`$, where $`w_i=P_i/\rho _i`$, and $`n_i=3(1+w_i)`$ (see, for example, Turner & White 1997; Peebles & Ratra 2003; or Daly & Guerra 2002). The deceleration parameter at the present epoch is $`q_o=(\ddot{a}a/\dot{a})_o=0.5\mathrm{\Omega }_i(1+3w_i)`$, when $`w_i`$ is time independent.

A direct, empirical determination of the acceleration of the universe as a function of redshift can be obtained from $`y(z)`$, without making any assumptions about the nature or evolution of the “dark energy.” This can be done using the equation (Daly 2002)

$$q(z)\ddot{a}a/\dot{a}^2=1+(1+z)(dy/dz)^1(d^2y/dz^2)$$
(6)
valid for $`k=0`$; if $`k0`$, another term $`[kr(1+z)/(1kr^2)](dr/dz)`$ must be added to the right hand side of equation (6).

Equation (6) depends only upon the Robertson-Walker line element and the relation $`(1+z)=a_o/a(t)`$. Thus, this expression for $`q(z)`$ is valid for any homogeneous, isotropic expanding universe in which $`(1+z)=a_o/a(t)`$, and is consequently quite general and can be compared with any model to account for the acceleration of the universe, as long as the model describes a homogeneous isotropic expanding universe with the standard relation between $`z`$ and $`a(t)`$. The Robertson-Walker line element is given by equation (4). A light ray emitted by a galaxy traveling to us along the radial coordinate $`r`$ has $`d\tau =d\theta =d\varphi =0`$. The increment is along the negative direction of $`dr`$ so eq. (4) with k=0 implies that $`a_odr=(1+z)dt`$, or $`(dz/dt)=a_o^1(1+z)(dr/dz)^1`$. Differentiating $`(1+z)=a_o/a(t)`$ with respect to time implies that $`\dot{a}=a_o(1+z)^2(dz/dt)`$. Substituting in for $`(dz/dt)`$, we find $`\dot{a}=(1+z)^1(dr/dz)^1`$. Differentiating again with respect to time, we find

$$\ddot{a}=(1+z)^2(dz/dt)(dr/dz)^1[1+(1+z)(dr/dz)^1(d^2r/dz^2)],$$
(7)
which simplifies to eq. (6) using the expressions given here, and the relation $`y(z)=H_0(a_or)`$. For $`k0`$, equation (4) implies that $`a_odr=(1+z)\sqrt{(1kr^2)}dt`$, or $`(dz/dt)=a_o^1(1+z)\sqrt{(1kr^2)}(dr/dz)^1`$. Differentiating $`(1+z)=a_o/a(t)`$ with respect to time implies that $`\dot{a}=a_o(1+z)^2(dz/dt)`$. Substituting in for $`(dz/dt)`$, we find $`\dot{a}=(1+z)^1\sqrt{(1kr^2)}(dr/dz)^1`$. Differentiating again with respect to time, we find $`\ddot{a}=(a_o)^1(1+z)^1(1kr^2)(dr/dz)^2[1+(1+z)(dr/dz)^1(d^2r/dz^2)]+(a_o)^1kr(dr/dz)^1,`$ thus $`q(z)=\ddot{a}a/\dot{a}^2=[1+(1+z)(dr/dz)^1(d^2r/dz^2)]+kr(1+z)(1kr^2)^1(dr/dz)`$.

Equation (6) can in principle be used to empirically determine the redshift at which the universe transitions from acceleration to deceleration without requiring assumptions regarding the nature and redshift evolution of the ”dark energy.” The supernova and radio galaxy data allow a determination of the dimensionless coordinate distance $`y`$ to each source, at redshift $`z`$. These data can then be used to determine $`dy/dz`$, and $`d^2y/dz^2`$; these can then be substituted into eq. (6) to determine q(z).

Since eqs. (5) and (6) are obtained without any assumptions regarding the mass-energy components of the universe or their redshift evolution, they can be used to directly determine the dimensionless acceleration rate $`E(z)=(dy/dz)^1`$, which contains important information on the “dark energy” and its redshift evolution, and to determine the dimensionless acceleration parameter $`q(z)`$ directly from measurements of $`y(z)`$.

In the determinations of $`y(z)`$ a value of $`_{}`$ must be adopted for the supernovae (see eq. 3), and a value of $`\kappa `$ must be adopted for radio galaxies (see eq. 1). There are not determined as a normalization using only low-redshift sources. They are determined by fitting all of the data and solving for the best fit values of these parameters. Fits to the supernovae data, the radio galaxy data, and the joint data set were run for a variety of cases (see Table 2), such as a universe with quintessence (Q), or a rolling scalar field (SF). There are enough supernovae that the value of $`_{}`$ changes very little for different fits to the supernovae data, and they change very little if the supernovae data are considered separately or in conjunction with the radio galaxy data. Thus, values of $`y`$ for supernovae do not change with the data set or model considered. The radio galaxy data best fit parameters for $`\kappa `$ and $`\beta `$ change slightly depending upon whether just the radio galaxies are considered, or whether the full data set of radio galaxies plus supernovae are included. Values of $`y_s`$ obtained for the best fit value of $`\kappa `$ using radio galaxies alone in a universe with quintessence are listed, as well as the values $`y_j`$ obtained using the best fit values of $`\kappa `$ and $`\beta `$ for fits to the full data set of radio galaxies and supernovae. These values are listed in Table 1 and are considered and compared in the analyses of $`E(z)`$ and $`q(z)`$.

## 4 The Numerical Differentiation Technique

The key problem in this approach, of course, is that it requires a numerical differentiation of typically noisy data, which is a cardinal sin for any empirical scientist. This, after all, is the reason why all standard cosmological tests (e.g., the Hubble diagram) consist of integrating the model equations to compare them with the measurements. An additional problem is posed by the sparse and/or uneven coverage of the redshift range(s) of interest. While a numerical differentiation of noisy data is in general not advisable, it is certainly possible, and if done properly (in a statistical sense), it can produce meaningful results within the limits of the available data.

Most numerical differentiation techniques explicitely or implicitely assume that the data can be locally represented by some smooth (differentiable) function, whose derivative is then defined analytically. Typically this local approximation is a low-order polynomial. Thus, estimation of derivatives is coupled with the estimation of the function representing the data themselves, in a self-consistent way. Measurement errors can then be propagated in the standard manner, leading to estimated uncertainties of the fitted function values as well as the derivatives. In our case, the function to be approximated, along with its first and second derivatives, is the dimensionless coordinate distance as a function of redshift, $`y(z)`$. The situation is simplified by the fact that the errors in $`z`$ are negligible in comparison to the errors in $`y`$, and thus the ordinary least-squares approach can be used.

There are three sources of errors when evaluating any function fits to noisy, finite data sets. First, the errors of the individual data points: the least-squares approach deals with them in a statistically optimal fasion, provided that the quoted error bars are truly representative, and that the deviations from the “true” underlying trend are drawn from a normal distribution. Second, if the fitted function is not a good approximation to the true trend, the results may be systematically biased. Locally, any function can be approximated as a polynomial (or as a Taylor series), and this becomes an issue of a having a sufficiently high fit order to account for the shape (the curvature) of the observed trend in the fitting interval. Finally, in any finite data set there will be some sample variance, i.e., a different draw of the same number of measurements from the same underlying trend, with the same errors, will produce slightly different results. The effects of the sample variance are minimised by having larger number of data points, and can be estimated numerically for any given sample.

We choose a simple powers-of-$`z`$ polynomial approach, in order to be maximally model-independent. In principle, other basis functions could be used, but we do not see any advantages of such an approach in a situation where the fits would be dominated by the noise and sparse sampling of the data. We always fit to $`y(z)`$, and then derive the first and second derivatives from the fit coefficients, and the local values of $`E(z)`$ and $`q(z)`$ using eqs. (5) and (6). Uncertainties of the fit coefficients are then propagated to derive the uncertainties in the fit values of $`y(z)`$, $`E(z)`$ and $`q(z)`$. The fit values are always evaluated on a redshift grid equally and densely spaced in either $`z`$ or $`\mathrm{log}z`$; this is just a matter of convenience, as the values and the quality of the fits are not affected.

A conceptually simplest approach would be to fit a polynomial to the entire data set. Unfortunately, low-order polynomials lack the flexibility to represent the actual shapes of underlying cosmological models, leading to seriously biased values of $`E(z)`$ and $`q(z)`$. The fits are (by design) optimised to fit the function $`(y)`$, and its derivatives are not constrained directly. Using higher order polynomials helps in recovering the mean shapes of these functions, but at the expense of greatly increased uncertainties, and typically with some oscillatory behavior, characteristic of high-order polynomial fits. For example, the $`q(z)`$ is generally a non-linear function of $`z`$, so the fits of an order $`>3`$ are needed; but in some cases, e.g., $`\mathrm{\Omega }_m=1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ cosmology, $`q(z)=const.`$, which higher order polynomials cannot reproduce very easily in shape, regardless of the increased errors for high-order terms.

A better method, which we adopted, is to fit the values of $`y(z)`$ locally, in some limited redshift window of $`\pm \mathrm{\Delta }z`$; within that interval, data points are fitted with the weights inversely proportional to the squares of their error bars. In addition, at each end of the fitting window we attach a Gaussian tapered region with a $`\sigma (z)=0.02`$, extending out to $`2\sigma `$; the enclosed data points in the tapered region have the weights lowered by the value of the Gaussian wing at that point. The purpose of this taper is to avoid fluctuations caused by individual data points entering and leaving the fitting window, at the expense of a slight increase in the resulting fit uncertainties (since the tapered points effectively get larger error bars). We established that the overall properties of the fits did not change. Finally, we require that there are at least 10 data points in each fit, and increment the window slightly if necessary.

The tradeoff in this technique is that larger fitting windows lead to more robust fits, at the expense of resolution in redshift and the introduction of the same problems which plague the global polynomial fits, as described above; while smaller fitting windows produce noisier fits because of a smaller number of enclosed data points. After some experimentation, we concluded that windows with $`\mathrm{\Delta }z0.4`$ seem to offer the optimal compromise, but we also perform fits with other window sizes.

After some experimentation, we decided to fit second order polynomials in each fitting window, as the minimal-assumption functions with defined second derivatives (needed to evaluate the $`q(z)`$), which can also accomodate any curvature in the data. We verified that using linear fits to obtain $`y^{}(z)`$ and thus $`E(z)`$ does not produce improved results, and that increasing the local fit order to 3 increases the formal errors without any significant benefits in terms of the fit quality and accuracy.

Specific details of the fitting procedure are as follows. Let the input data be $`(z_i,y_i,\mathrm{\Delta }y_i)`$. The fitting weights are computed in a standard fashion as $`1/\mathrm{\Delta }y_i^2`$. The fits are evaluated on an output reshift grid $`z_j`$, typically ranging from 0.01 to 1.7, with a spacing of 0.01 or 0.005. This is simply a choice of convenience, since the fits can be evaluated anywhere in the redshift regime covered by the data, and there is no reason to do it, say, just at the values of the input $`z_i`$. We note that since our fitting windows $`\mathrm{\Delta }z`$ are generally much larger than the output grid spacing, the adjacent output fit values are $`not`$ independent. For each output point, $`z_j=z_0`$, we select the input data in the corresponding fitting window as described above, ranging from some $`i=i_{min}`$ to some $`i=i_{max}`$; the fitting weights for points in the Gaussian taper regions are adjusted appropriately. We perform the “centered” fits by introducing the independent variable $`x_i=z_iz_0`$, and fit second order polynomials to $`y(x)`$ in the range from $`i_{min}`$ to $`i_{max}`$. We use the routine $`fit`$ from Numerical Recipes (Press et al. 1992). The output are the fitting coefficients $`A,B,C`$, where $`y(x)=A+Bx+Cx^2`$. Then at the centered fit value of $`x=0`$, the best fit value for $`y(z)=A`$, the best fit value of its first derivative is $`dy/dzdy/dx=B`$, and the best fit value of the second derivative is $`d^2y/dz^2d^2y/dx^2=2C`$. Thus, from eq. (5), $`E(z)=1/B`$, and from eq. (6), $`q(z)=1(1+z)\times 2C/B`$. The routine $`fit`$ also returns the covariance matrix, whose diagonal elements give the uncertainties of the coefficients $`A,B,C`$, and their uncertainties are easily propagated to the fitting uncertainties of $`E(z)`$ and $`q(z)`$.

We emphasize that our goal is not to evaluate a number of independent measurements of $`E(z)`$ and $`q(z)`$ in the redshift range of interest (the data in hand are not sufficient to do that for an interesting number of points), but instead to outline the global trends presented by the data. Thus, we use relatively broad fitting windows with our sliding window fit methodology. The price we pay is the strong correlation of fitted values on our densely spaced output grid, and these should not be taken as independent measurements, but really as outlines of the global trends. In our model-independent approach there is an implicit (and reasonable) assumption that the function $`y(z)`$ is smoothly changing, with only a modest local curvature. This implies that there is a useful “nonlocal” information present in the data, which is captured by our extended fitting windows, at a price of having only 2 or 3 “independent” measurements across our full redshift range. With a richer input data set we could make the sliding windows smaller, and increase the number of “independent” measurements of $`E(z)`$ and $`q(z)`$ (i.e., our redshift resolution of the measured trends).

An alternative approach would be to bin the data in a modest number of redshift bins, and perform an independent fit in each bin. An advantage of such a technique would be a set of truly independent measurements of $`E(z)`$ and $`q(z)`$, but at a price of an increased noise, since the information present by points outside each bin boundary would be lost. This would also lead to implied discontinuities in the fit values of $`y(z)`$ at the bin edges, leading to implied infinite derivatives, which is obviously unphysical. Thus, we opt for a sliding window fit approach, which is essentially a flexible and robust smoothing technique, with a clear caveat that the output fit values are not independent within the used $`\mathrm{\Delta }z`$ range.

In exploring the numerical fitting and differentiation methodology, we used artificial data sets with known, built-in cosmologies, in order to evaluate the accuracy of the derived fits for $`E(z)`$ and $`q(z)`$. For most part, we generated artificial data sets mimicking what is expected from SN measurements by the SNAP satellite (see, e.g., Aldering et al. 2003, or http://snap.lbl.gov/), namely a set of 2000 measurements in the redshift interval from 0.1 to 1.7, with combined (measurement + intrinsic) scatter of 7% in dimensionless coordinate distances. The redshift distribution function was taken to be proportional to the volume element divided by the redshift, which roughly represents a combination of the expected SN rate history and the SNAP selection function. For most tests, we assumed the standard Friedmann-Lemaitre cosmology with $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. For some of the tests we changed some of these assumptions (the number of the data points, the relative errors, or the underlying cosmology).

The overall process is illustrated in Fig. 3, on an example of our pseudo-SNAP data set. The global trends for $`y(z)`$, $`E(z)`$, and $`q(z)`$ are reproduced well, with the bias (systematic offsets) comparable to the fit errors. The errors increase going towards the higher derivatives and near the edges of the redshift intervals, as may be expected.

The effects of different fitting windows on the derived values $`y(z)`$ and $`E(z)`$ are illustrated in Fig. 4. As expected, smaller values of $`\mathrm{\Delta }z`$ lead to noisier fits, but the overall trends agree within the errors. For the present RG+SN data set, we use $`\mathrm{\Delta }z`$ of 0.4 or 0.6 in what follows.

Finally, we address the issue of the sample variance. For our pseudo-SNAP data, we simply generate a number of different realizations of the data set, using the same assumptions. In order to estimate the effects for our RG+SN data set, we do the following. We assume an underlying cosmology, viz., $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. Then, we replace each $`y(z)`$ measurement with the value for this cosmology, perturbed by a random amount drawn from a Gaussian distribution with the $`\sigma (y)`$ equal to the quoted error bar. We generate a number of such pseudo-RGSN data sets, and perform the fits on each. The results are shown in Fig. 5. For the RG+SN data set, the sample variance effects are comparable to or smaller than the measurement errors for $`y(z)`$; and comparable to the fitting uncertainties (which derive from the random errors of the data) for $`E(z)`$. As expected, the sample variance effects for the pseudo-SNAP data set, which has many more data points, are effectively negligible.

Estimates of the sample variance errors are made only rarely in the published literature. Our tests suggest that for many real-life data sets in cosmology, these errors can be easily comparable to the fitting uncertainties which derive from the random measurement errors, and thus in many cases the quoted confidence intervals may be underestimating the total uncertainties.

## 5 The Initial Results for E(z) and q(z)

Using the procedure described above, the function $`E(z)`$ was obtained for the full RG+SN data set using eq. (5), and is shown in Fig. 6. We used the values of $`y_j`$ for RGs listed in Table 1, and the values for $`y`$ for SNe listed in Tables 3 and 4; a total of 78 SNe were used including the average values of $`y`$ for the 14 SNe listed in Table 4, and values of $`y`$ for the remaining 64 SNe listed in Table 3. The results are remarkably close to the currently popular “concordance” Friedmann-Lemaitre model with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`\mathrm{\Omega }_m=0.3`$. We note, however, that in our analysis we did not assume that the universe is described by a Friedmann-Lemaitre model at all.

As an internal consistency test, we computed the fits for the RG and SN samples separately, using the values $`y_s`$ for the radio galaxies. The results are shown in Fig. 7. It is clear that the values of $`y`$ obtained for radio galaxies and supernovae agree very well for the redshifts where the data sets overlap. And, in the redshift range where the two samples overlap, the independent determinations of $`E(z)`$ agree to within 1-$`\sigma `$ (joint errors) or better. It is notable that the RG data alone are consistent with a constant $`E(z)1.1`$ for $`z0.41.8`$, although the error bars are large. If this trend remains as the error bars decrease with better and more extensive data sets, it could be indicative of an actual cosmological trend, or (perhaps more likely) some evolutionary effect or bias in the data. This can not be sorted out until more data are avialable.

Finally, we show in Fig. 8 what is probably the first direct estimate of $`q(z)`$, obtained using eq. (6). The data, folded through our analysis procedure, are fully consistent with the “concordance” model with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`\mathrm{\Omega }_m=0.3`$, and suggest that indeed the universe transitions from decceleration to acceleration at $`z_T\text{ }>\text{ }0.3`$, with a best fit value of $`z_T0.45`$. Again, we note that no assumptions about the cosmological model have been made in deriving this trend. This is a preliminary result, and we are clearly limited by the available data at this time. Our purpose here is mainly to illustrate the method, but even so, the results are very encouraging.

## 6 Implications of E(z) and q(z) for the Properties of the Dark Energy

The acceleration parameter is $`q(z)=\ddot{a}a/\dot{a}^2=(\ddot{a}/a)(\dot{a}/a)^2`$. The acceleration of the universe is described by

$$(\ddot{a}/a)=\frac{4\pi G}{3}(\rho _i+3p_i)=\frac{4\pi G}{3}\rho _i(1+3w_i)$$
(8)
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`$. In addition, for $`k=0`$,

$$(\dot{a}/a)^2=\frac{8\pi G}{3}\rho _i.$$
(9)
Thus, $`q(z)=0.5\rho _i(1+3w_i)/\rho _i`$ when $`k=0`$, and $`E^2(z)=\rho _i/\rho _{co}`$, where $`\rho _{co}`$ is the critical density at the present epoch, $`\rho _{co}(3/8\pi G)H_o^2`$ = $`\rho _{mo}+\rho _{Eo}`$, and the present epoch mean mass-energy density of non-relativistic matter and dark energy are $`\rho _{mo}`$ and $`\rho _{Eo}`$ respectively. Non-relativistic matter evolves with redshift as $`(1+z)^3`$. Let the dark energy evolve with redshift as $`f_E(z)`$, then, it is easy to show that $`E^2(z)=\mathrm{\Omega }_m(1+z)^3+(1\mathrm{\Omega }_m)f_E(z)`$, where $`\mathrm{\Omega }_m=\rho _{mo}/\rho _{co}`$, and $`\mathrm{\Omega }_E=\rho _{Eo}/\rho _{co}=1\mathrm{\Omega }_m`$. Hence, if the current contribution of non-relativistic matter $`\mathrm{\Omega }_m`$ can be determined, then $`E(z)`$ can be used to determine the redshift evolution of the dark energy $`f_E(z)`$. Similarly, $`q(z)=0.5[\mathrm{\Omega }_m+(1\mathrm{\Omega }_m)(1+3w)(1+z)^3f_E(z)]/[\mathrm{\Omega }_m+(1\mathrm{\Omega }_m)(1+z)^3f_E(z)]`$.

Now, for quintessence, $`\rho _i=\rho _{i,o}(1+z)_i^n`$ when $`w_i`$ is constant, where $`n_i=3(1+w_i)`$. This follows from the mass-energy conservation of each component, which implies

$$\dot{\rho _i}=3(\rho _i+p_i)(\dot{a}/a)$$
(10)
(e.g. Peebles 1993). When the equation of state $`w_i`$ does not change with time, 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)`$.

With two important components at low redshift, a non-relativistic component $`\rho _m`$ that includes baryons and the dark matter in galaxies and clusters of galaxies, and a dark energy component $`\rho _E`$, we have $`\rho _m=\rho _{mo}(1+z)^3`$ and $`\rho _E=\rho _{Eo}(1+z)^n`$. Now, at zero redshift the total density is equal to the critical density $`\rho _{co}=\rho _{mo}+\rho _{Eo}`$ and since $`\rho _{mo}/\rho _{oc}\mathrm{\Omega }_m`$, the acceleration parameter may be written

$$q(z)=0.5\left(\frac{\mathrm{\Omega }_m+(1\mathrm{\Omega }_m)(1+3w)(1+z)^{3w}}{\mathrm{\Omega }_m+(1\mathrm{\Omega }_m)(1+z)^{3w}}\right).$$
(11)
The universe is decelerating when the sign of $`q(z)`$ is positive. The sign of the denominator is always positive, and the numerator may be either positive or negative depending upon the value of $`w`$. The universe will go from a state of acceleration to a state of deceleration if the dark energy has properties like that of quintessence when the numerator of equation (8) is positive, which occurs a the transition redshift $`z_T`$ given by

$$z_T=\left[\frac{\mathrm{\Omega }_m}{(1+3w)(1\mathrm{\Omega }_m)}\right]^{1/(3w)}1.$$
(12)
This transition redshift is plotted as a function of the equation of state $`w`$ in Figure 9. Clearly, as the transition redshift increases, the value of $`\mathrm{\Omega }_m`$ must decrease, or the equation of state $`w`$ exhibit redshift evolution.

Similarly, for quintessence, $`E^2(z)=(1+z)^3[\mathrm{\Omega }_m+(1\mathrm{\Omega }_m)(1+z)^{3w}]`$. Some lines representing quintessence with $`w=1`$ (i.e. a cosmological constant) are included in the figures.

Our preliminary results give a limit to the transition redshift $`z_T\text{ }>\text{ }0.3`$, with the best fit estimate $`z_T0.45`$. Assuming $`\mathrm{\Omega }_m=0.3`$, these translate to $`w\text{ }<\text{ }0.55`$ and $`2.3\text{ }<\text{ }w\text{ }<\text{ }0.65`$ (see Figure 9). With better data sets in the future, we should be able to improve on these limits.

Note that $`q(z)`$ is an important input into the Statefinder diagnostic presented by Sahni et al. (2002) and Alam et al. (2003); these authors discuss one way in which $`q(z)`$ may be used to determine $`w(z)`$ and the Statefinder pair (r,s).

## 7 Summary and Discussion

We presented here a set of dimensionless coordinate distances for 20 RGs, spanning the redshift range 0.43 to 1.79 (with one source at $`z=0.056`$). These measurements supplement and extend to the cosmologically interesting redshift range the distances available for SNe, which currently reach only to $`z=0.97`$ (with one source at $`z=1.70`$).

The determination of the dimensionless coordinate distances to RGs and SNe are completely independent, and are based on completely different physics. Yet, the two data sets agree very well in the overlap redshift range, as shown here, and as shown by previous work (e.g., Podariu et al. 2002; Daly & Guerra 2002; Guerra, Daly, & Wan 2000). This is very encouraging: there is a great value in being able to measure the same physical quantity (here the coordinate distances as a function of redshift) using different and independent tracers and methods. The general agreement we see between the SN and RG data sets suggests that neither method is dominated by some substantial, as yet unknown systematic errors. Together, the two data sets can be used in cosmological tests with a greater power than each data set separately.

The dimensionless coordinate distances $`y(z)`$ can be used to empirically determine the dimensionless expansion and deceleration rates as functions of redshift, $`E(z)`$ and $`q(z)`$, without assuming any particular cosmological model. While the traditional cosmological tests integrate the expressions for these functions provided by the models (e.g., the standard Friedmann-Lemaitre models) and determine the model parameters from the fits, we develop a complementary procedure whereby these functions can be derived directly from the data by differentiating the $`y(z)`$ trend. We apply a particular, simple numerical procedure to this task, and derive the trends of both $`E(z)`$ and – for the first time – $`q(z)`$ directly from the data.

Our estimates of $`E(z)`$ are in an excellent agreement with those obtained from other methods, e.g., the CMBR fluctuations,large-scale structure, high-$`z`$ SNe and radio galaxies using traditional analyses, etc., even through they are obtained in a completely different and independent manner. In particular, the data are consistent with the “concordance” cosmology, i.e., Friedmann-Lemaitre models with $`\mathrm{\Omega }_m=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. While these results are clearly very preliminary, and meant primarily to illustrate the method, the good agreement with other approaches is very encouraging.

We are currently limited by the amount and quality of the available data for both SNe and RGs. Nevertheless, there are great prospects for advances in precision cosmology, e.g., large sets of high-quality measurements of SNe from the SNAP satellite (e.g., Aldering et al. 2003), or from large ground-based experiments such as the ESSENCE (Stubbs 2002; see also http://www.ctio.noao.edu/wproject/) or LSST in the future (Tyson et al. 2002, 2003). Such data sets could certainly support differentiation of distance vs. redshift trends, leading to considerably more robust direct determinations of the expansion and acceleration rates as functions of redshift. In addition, new radio galaxy data is being obtained.

As the observational situation improves, direct estimates of $`E(z)`$ and $`q(z)`$ can be used to understand the properties and redshift evolution of different flavors of dark energy, determine the redshift at which the universe transitions from acceleration to deceleration, and may help elucidate any systematic errors that might be lurking in the RG or SN methods of constraining cosmological parameters.

It is a pleasure to thank Megan Donahue, Eddie Guerra, Matt Mory, Chris O’Dea, Paddy Padmanabhan, Bharat Ratra, and Varun Sahni for helpful comments and discussions, Saul Perlmutter for sending electronic files of supernova data, and Adam Reiss for providing us with supernova data modified to be on the same scale as that provided by Saul Perlmutter. This work was supported in part by the U. S. National Science Foundation under grants AST-0096077 and AST-0206002, and Penn State University (RAD), and by the Ajax Foundation (SGD).


text:
# Discovery of a New Pulsar Wind Nebula in the Large Magellanic Cloud

## 1. Introduction

The Magellanic Clouds are ideal regions for studying supernova remnants (SNRs), the known distance and low extinction obviating many of the frustrations associated with studying Galactic SNRs. While this has proved useful in studying the shell emission produced by the interaction of SNRs with the ambient medium, just two Magellanic SNRs, B0540–693 (Manchester, Staveley-Smith, & Kesteven (1993)) and N157B (Wang et al. (2001)), are known to harbor a young energetic pulsar and its surrounding pulsar wind nebula (PWN). The lack of other PWNe in the Magellanic Clouds is likely due to a lack of spatial resolution — such resolution is only now being applied to Magellanic SNRs in both the radio and X-ray bands, and so the prospects are good for identifying new “composite” SNRs, containing both outer shells and central PWNe. Such sources are of particular interest, as the shell and PWN together provide unique constraints on the properties of the system.

We here present data which demonstrate that B0453–685 in the LMC is a composite SNR. This Letter primarily discusses the central PWN; interpretation of the surrounding shell will be the focus of a subsequent paper (Hendrick et al. (2003)).

## 2. Observations and Results

Radio observations of B0453–685 were carried out with the Australia Telescope Compact Array (ATCA). Observations were made simultaneously at 1.4 and 2.4 GHz with a bandwidth of 128 MHz at each frequency, using the 1.5G1111.5G is a non-standard ATCA configuration, and has successive spacings between the six antennas of 30.6, 535.7, 306.1, 413.3 and 3015.3 meters. and 6.0C configurations on 2002 Jul 23 and 2002 Aug 26 respectively. A 38.6-ks X-ray observation of the same field was carried out using the ACIS-S3 detector of the Chandra X-ray Observatory on 2001 Dec 18.

For each radio observation, the time on source was $``$11 hr. Antenna gains and polarization leakages were calibrated using observations of PKS B0252–712, while the flux density scale was tied to PKS B1934–638. The radio data were edited, calibrated and imaged using standard techniques;222See http://www.atnf.csiro.au/computing/software/miriad . the resulting resolution at each frequency is listed in Table 1.333The 1.4-GHz image was formed using all antennas, while the 2.4-GHz image excluded the 6th antenna so as to not overresolve the source. Analysis of the Chandra data will be reported elsewhere (Hendrick et al. (2003)).

### 2.1. Imaging and Polarimetry

Radio and X-ray images of SNR B0453–685 are shown in Figure 1. At radio frequencies, the source is clearly resolved into two components: a bright central core, embedded in an approximately circular limb-brightened shell. Table 1 lists the flux densities of these two components at 1.4 and 2.4 GHz after applying a correction for the local background. Fitting an ellipse to the perimeter of the shell component, we find that if one excludes the bright protrusion to the southwest, the shell can be fitted by a circle of diameter $`120^{\prime \prime }\pm 4^{\prime \prime }`$, centered at RA (J2000) $`04^\mathrm{h}53^\mathrm{m}37\stackrel{\mathrm{s}}{\mathrm{.}}2`$, Decl. (J2000) $`68^{}29^{}30^{\prime \prime }`$ (with an uncertainty $`\pm 2^{\prime \prime }`$). For a distance $`50d_0`$ kpc, the diameter of the shell is $`(29\pm 1)d_0`$ pc.

Figure 2 shows that the central radio core is elongated along a position angle $`45^{}`$, with dimensions $`30^{\prime \prime }\times 20^{\prime \prime }`$ (i.e. a spatial size $`7.3d_0`$ pc $`\times `$ $`4.8d_0`$ pc). The core is brightest in the NE, where it is dominated by a slightly extended region, with peak position (J2000) RA $`04^\mathrm{h}53^\mathrm{m}38\stackrel{\mathrm{s}}{\mathrm{.}}6`$, Dec. $`68^{}29^{}21\stackrel{}{\mathrm{.}}7`$ ($`\pm 0\stackrel{}{\mathrm{.}}5`$ in each coordinate). The core shows significant linear polarization, of mean fractional intensity 6% at 1.4 GHz and 8% at 2.4 GHz, as shown in Figure 2. The brightest region of the core is offset from the center of the shell by $`11^{\prime \prime }\pm 2^{\prime \prime }`$. To the SW, the surface brightness of the core fades with increasing distance from the peak.

The third panel of Figure 1 demonstrates that the morphology of B0453–685 in the soft X-ray band shows a broadly similar structure to that seen in the radio, again revealing a central core within a limb-brightened shell. In the hard X-ray band (lower panel of Fig 1), only the central core is visible. The extent of this region is $`14^{\prime \prime }\times 7^{\prime \prime }`$ ($`3.4d_0`$ pc $`\times `$ $`1.7d_0`$ pc), less than half that seen in the radio. Figure 2 demonstrates that the positions of the peak emission in the X-ray and radio cores coincide to within the resolution of the observations.

### 2.2. Spectral Analysis

We have carried out a spectral index study of the radio data by smoothing the 1.4 and 2.4 GHz images to a common resolution of $`15^{\prime \prime }`$ and then applying spectral tomography to these maps (Katz-Stone & Rudnick (1997)). We find that the central core has a flatter radio spectrum than the surrounding shell: tomography implies $`\alpha =0.10\pm 0.05`$ for the core and $`\alpha =0.5\pm 0.1`$ for the shell ($`S_\nu \nu ^\alpha `$). Assuming $`\alpha =0.10`$ for the core between $`10^7`$ and $`10^{11}`$ Hz, the radio luminosity of this component is $`L_R2d_0^2\times 10^{34}`$ ergs s-1.

At X-ray energies, we have determined the spectral properties of the different components of the SNR by extracting events from appropriate source and background regions and fitting to the resulting spectra, as will be described in detail by Hendrick et al. (2003). To briefly summarize their results, the shell component is well-fitted by the absorbed thermal spectrum expected in the Sedov phase of SNR evolution, while the central core is best fit by an absorbed power law.444 Bremsstrahlung ($`kT4`$ keV) can also fit the central core, but the absence of lines demands either absurdly low abundances or an extremely underionized plasma (ionization timescale $`n_et<100`$ yr cm-3). The best-fit parameters for the core are a foreground absorbing column $`N_H=(1.3\pm 0.2)\times 10^{21}`$ cm-2, a spectral index $`\alpha =0.9\pm 0.4`$ and an unabsorbed flux density at 1 keV $`f_X=(6\pm 1)\times 10^{14}`$ erg cm-2 s-1 keV-1 (uncertainties quoted at 90% confidence). The corresponding X-ray luminosity over the energy range 0.5–10 keV is $`L_X6d_0^2\times 10^{34}`$ erg s-1. The spectral fits for the surrounding shell imply an age $`\tau 13`$ kyr, a swept-up mass $`M_{sw}280M_{}`$, an upstream density $`n_00.4`$ cm-3 and an initial explosion energy $`E_05\times 10^{50}`$ ergs (Hendrick et al. (2003)).

## 3. Discussion

The core embedded in B0453–685 is centrally located, polarized, has a filled-center morphology and a flat radio spectrum, and is coincident with a smaller non-thermal X-ray nebula. These properties demonstrate that B0453–685 is a composite SNR, containing a synchrotron-emitting PWN powered by a neutron star embedded in a shell SNR.

### 3.1. Properties of the Central Pulsar

Within the core, no radio or X-ray point source corresponding to the central pulsar is apparent. We can calculate upper limits on the flux from such a source as follows. At radio wavelengths, we made images of the field using only the longest ($`>3000`$ meters) five baselines in each configuration. At both 1.4 and 2.4 GHz, the brightest region of the core was still detected in these images. At the peak of nebular emission, the corresponding upper limit on a point source is 3 mJy and 0.4 mJy at 1.4 and 2.4 GHz, respectively.

In the X-ray band, we employed the SHERPA package to carry out a 2D fit to the nebular morphology. We fit the core with three components: an elliptical gaussian to define the overall shape, an unresolved source to represent an embedded pulsar, and a constant offset to account for background. Using these three components, we find that in the energy range 0.3–10 keV, the best-fit model contains 58 counts in an unresolved central source. While the apparent significance of this source seems high, we obtained equally good or better fits by replacing the unresolved component by a resolved Gaussian. This suggests a complex PWN morphology and only marginal evidence for a point source. Therefore, in subsequent discussion we adopt 58 counts as an upper limit, corresponding to a count-rate $`<1.5\times 10^3`$ cts sec-1. Assuming a column density $`N_H=1.3\times 10^{21}`$ cm-2 and a power law spectrum with a photon index $`\mathrm{\Gamma }=1.5`$ (typical for a young pulsar), this corresponds to an unabsorbed X-ray luminosity (0.5–10 keV) $`L_X<6d_0^2\times 10^{33}`$ ergs s-1.

These radio and X-ray upper limits on emission from a central pulsar are not especially constraining. Our radio limit is $`40`$ times above the upper limit for radio pulsations at this position from the survey of Manchester et al. (2003). Several young pulsars have been recently identified with a point-source X-ray luminosity $`L_X10^{32}10^{33}`$ ergs s-1 (Murray et al. (2002); Hughes et al. (2003)), well below the limit seen here.

In the absence of a direct detection, a first guess as to the properties of the pulsar in B0453–685 can come from comparison with other systems. There are $`20`$ other shell SNRs known to contain a central X-ray/radio PWN; of these other sources, the properties of B0453–685 most closely resemble those of the Vela SNR and of G0.9+0.1, as demonstrated in Table 2. The similarity of these three systems argues that SNR B0453–685 is most likely powered by a “Vela-like” pulsar (e.g. Kramer et al. (2003)), with a spin-period $`P100`$ ms, a surface magnetic field $`B3\times 10^{12}`$ G and a spin-down luminosity $`\dot{E}10^{37}`$ ergs s-1.

In cases in which a PWN has no identified pulsar, considerable further effort has been invested to infer the pulsar’s properties (i.e. current spin period $`P`$, initial period $`P_0`$, period derivative $`\dot{P}`$ and surface magnetic field $`B`$) from those of its PWN. The age, spin-down luminosity and surface magnetic field of the system are, respectively:

$$\tau =\frac{P}{(n1)\dot{P}}\left[1\left(\frac{P_0}{P}\right)^{n1}\right],\dot{E}=4\pi ^2I\frac{\dot{P}}{P^3},B=3.2\times 10^{19}(P\dot{P})^{1/2}\mathrm{Gauss},$$
(1)
where $`I10^{45}`$ g cm2 is the star’s moment of inertia, $`n2P\ddot{P}/\dot{P}^2`$ is the “braking index”, and $`P`$ is in seconds. Typically one assumes $`n=3`$, $`P_0=0`$ and $`L_X=\eta \dot{E}`$ (where $`\eta `$ is an assumed efficiency factor). $`L_X`$ and $`\tau `$ are determined from observations, from which estimates for $`P`$, $`\dot{P}`$ and $`B`$ can then be derived using Equations (1).

Here we propose an alternative approach, where rather than assuming an (unphysical) value of $`P_0`$, we use the above expressions to eliminate $`P`$ and $`\dot{P}`$, and thus obtain a relation between $`B`$ and $`P_0`$. Assuming that $`B`$ is constant with time, the resulting function has absolute maxima in both $`B`$ and $`P_0`$. The robust upper limits on the surface magnetic field and the initial spin period of the neutron star which result are:

$$B<3.2\times 10^{19}\frac{U}{\tau (n1)}\mathrm{Gauss};P_0<U(\tau [n+1])^{1/2}\left(\frac{2}{n+1}\right)^{\frac{1}{n1}},$$
(2)
where $`U^2=P^3/\dot{P}=4\pi ^2I\eta /L_X`$. These limits are reasonably insensitive to the value of $`n`$ assumed; for $`2<n<3`$ (encompassing the range for the four pulsars for which $`n`$ has been directly measured; e.g. Cusumano, Massaro, & Mineo (2003)), we find $`B<(3.16.3)\eta ^{1/2}\times 10^{13}`$ G and $`P_0<(448488)\eta ^{1/2}`$ ms for this system. For other PWNe, a solid upper limit is $`\eta <0.05`$; we have argued above that the pulsar here is “Vela-like”, for which $`\eta 0.010.001`$ is more reasonable (Possenti et al. (2002)). We thus conclude that this pulsar has a surface magnetic field $`B<6\times 10^{12}`$ G and an initial spin period $`P_0<50`$ ms. These values are consistent with “typical” radio pulsars such as the Crab ($`B=4\times 10^{12}`$ G, $`P_0=19`$ ms) and Vela ($`B=3\times 10^{12}`$ G) pulsars, but not with the emerging class of young pulsars which are highly magnetized and/or slow initial rotators, such as PSRs J1846–0258 ($`B=5\times 10^{13}`$ G; Gotthelf et al. (2000)) or J1210–5226 ($`P_0400`$ ms; Pavlov et al. (2002)).

### 3.2. Evolutionary State of the PWN

Chevalier (1998) divides PWN evolution into successive phases. The PWN first expands supersonically in the SNR interior. The PWN then collides with the SNR reverse shock; this compresses the PWN and causes it to expand more slowly. The pulsar’s motion can later become supersonic in the shocked SNR interior, and it then drives a bow-shock.

If we assume that the pulsar spin-down luminosity is constant in the initial supersonically expanding phase, the radius of the PWN evolves such that $`R_{PWN}=0.839\times (\dot{E}\tau /E_0)^{1/5}V_0\tau `$, where $`V_0`$ was the initial expansion velocity of the SNR (van der Swaluw et al. (2001)). Here we observe $`R_{PWN}3`$ pc, $`\tau 13`$ kyr and $`E_05\times 10^{50}`$ ergs, and we adopt $`\dot{E}=E_{37}\times 10^{37}`$ ergs s-1. Thus if the PWN is still in this phase of evolution, we require $`V_0700E_{37}^{1/5}`$ km s-1, implying an enormous ejected mass $`M_{ej}=2E_0/V_0^2100E_{37}^{2/5}M_{}`$. We thus conclude that the PWN is almost certainly not expanding supersonically, having a radius much smaller than expected in this interpretation.

The brightest radio and X-ray emission from the PWN is offset from the SNR’s geometric center by $`l=(2.7\pm 0.5)d_0`$ pc. Since this peak most likely marks the current position of the unseen pulsar, a possible interpretation is that the pulsar is moving away from its birth site with a transverse velocity $`V_T=l/\tau =200\pm 40`$ km s-1. The radio PWN is clearly elongated, with a fading tail pointing back along the implied direction of motion. While this is suggestive of a bow shock morphology, this requires that the pulsar be moving supersonically through shocked ejecta. In the Sedov solution, this condition is met only when $`l/R_{SNR}0.7`$ (van der Swaluw, Achterberg, & Gallant (1998)). Unless the pulsar’s motion is at less than $`20^{}`$ to the line-of-sight (implying a true velocity $`V>700`$ km s-1), the pulsar in this system is still too close to the center of its SNR to meet this requirement.

The remaining possibility is that this PWN has undergone a reverse shock interaction with its SNR, and is now expanding subsonically. Such an interaction begins at an age $`10`$ kyr when $`M_{sw}10M_{ej}`$ (Chevalier (1998)), consistent with the properties of this source. The elongated morphology and comparatively small radius seen for the radio PWN are simply accounted for by reverse-shock compression; the same properties have been similarly interpreted in the other SNRs listed in Table 2 (Blondin, Chevalier, & Frierson (2001)).

## 4. Conclusions

New radio and X-ray observations of SNR B0453–685 in the LMC demonstrate this source to consist of both an outer shell and a central pulsar-powered nebula. We have used the properties of this SNR to infer that the unseen central engine is a typical young radio pulsar, with an initial period $`P_0<50`$ ms, a current period $`P100`$ ms, a surface field $`B<6\times 10^{12}`$ G and a spin-down luminosity $`\dot{E}10^{37}`$ erg s-1. The small and elongated radio nebula results from compression of the PWN by the SNR reverse shock.

This study demonstrates that several useful constraints on the properties of an unseen pulsar can be inferred from those of its associated PWN and SNR, provided that good estimates for the system’s distance and age are available. With many new studies of PWNe and composite SNRs now emerging, such an approach can be applied to many other sources.

We thank Bob Sault for allowing us to carry out the ATCA observations remotely from Harvard University, Dick Manchester for information on pulsar surveys of the LMC, and Eric van der Swaluw and Patrick Slane for useful discussions. The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. This work was supported by NASA through grants GO1-2075X and NAG5-13032.