text:
# 𝐾→𝜋𝜋 decay amplitudes from the lattice Presented by C. Kim. This work was supported by in part by the U. S. Department of Energy.

## 1 Introduction

The techniques of lattice gauge theory offer the possibility of accurate, non-perturbative calculation of the masses and matrix elements predicted by low-energy QCD. Of particular interest are the properties of the K-meson system and the low-energy, weak matrix elements that are required to understand Kaon decays, especially those which enter the $`K\pi \pi `$ amplitudes which violate CP symmetry. Presently well-developed methods permit the calculation of $`K`$ to vacuum and $`K`$ to $`\pi `$ weak matrix elements. When combined with chiral perturbation theory and the accurate chiral symmetry of the domain wall fermion formulation, lattice methods can be used to evaluate all the matrix elements that contribute to the weak $`K\pi \pi `$ amplitudes in the chiral limit. In this paper we describe a method which may provide a practical approach to the direct lattice calculation of the desired, on-shell, $`K\pi \pi `$ amplitudes without the use of chiral perturbation theory.

## 2 Finite volume decay amplitudes

We begin by reviewing the finite-volume approach developed by Lellouch and Lüscher. To compute the matrix element of an operator $`O`$ between two physical states $`|A`$ and $`|B`$, one typically evaluates the Euclidean Green’s function:

$$G(t_1,t_2,t_3)=\mathrm{Tr}\left\{B(t_1)O(t_2)A(t_3)e^{HT}\right\}$$

in the limit $`Tt_1t_2t_3`$. Here the operators $`A`$ and $`B`$ are chosen to carry quantum numbers which insure that the desired states, $`|A`$ and $`|B`$ are each the lowest energy states that can be created from the vacuum by $`A`$ and $`B`$, respectively. As demonstrated by Maiani and Testa, this method fails when one of the desired states is a two-particle state above threshold. In this case the desired final state cannot be uniquely selected by its quantum numbers and the Euclidean time limit will give an unwanted, threshold amplitude.

In a recent paper, Lellouch and Lüscher propose a method to circumvent this difficulty by tuning the finite lattice volume so that the first excited two-particle state allowed in that volume has the energy of the initial decaying particle. They work out an explicit formula which relates the resulting finite-volume matrix element to the desired infinite volume decay amplitude. The same result is obtained by different methods by Lin et al., offering further insight. As proposed, this method faces serious practical difficulties given our present computational capabilities. First, a large spatial volume of $`(6\mathrm{Fermi})^3`$ is required if the first-excited $`\pi \pi `$ state is to have an energy matching the Kaon mass. Second, the extraction of the next-leading exponential in the large-time limit is numerically difficult, especially if an accurate result is desired. Finally, while the matrix element of the $`I=2`$ final state can be obtained from the coefficient of such a next-leading exponential, the more important $`I=0`$ state also receives a vacuum contribution implying that the desired amplitude is the third term in an expansion in increasingly small exponentials.

We propose to overcome these difficulties by imposing an anti-periodic boundary condition on the pions, say in the $`z`$-direction, while using the usual periodic boundary conditions in the other two directions. With this choice the lowest energy state of a single pion has momentum $`p_z=\frac{\pi }{L}`$ instead of zero. With such a choice of boundary conditions, the lowest energy state with $`I=2`$ will be two pions, each with momentum $`p_z\pm \frac{\pi }{L}`$. Thus, the matrix element of such a state can be extracted from the leading exponential in the Euclidean Green’s function. Since the vacuum state is unaffected by these boundary conditions, our $`I=0`$ two-pion state will show the next-leading large-time behavior and a vacuum subtraction will be necessary for calculation of $`\mathrm{\Delta }I=\frac{1}{2}`$ amplitudes.

The condition which insures that the two-pion state has the relative momentum required for a physical $`K\pi \pi `$ decay, now requires $`L_z3`$ Fermi, only one-half the lattice extent needed for the periodic case. The other two spatial dimensions are not constrained by these kinematics and might be chosen to be similar $`L_x=L_y=L_z`$ giving a reduced box size of $`(3\mathrm{Fermi})^3`$. Thus, the proposed method decreases the required spatial size from $`(6\mathrm{Fermi})^3`$ to $`(3\mathrm{Fermi})^3`$ and promotes the amplitude of interest to one of leading ($`I=2`$) or next-leading ($`I=0`$) large-time behavior. Since only periodic and cubically symmetric periodic boundary conditions were discussed in the original finite volume analysis on which the Lellouch-Lüscher paper is based, we must recompute Lüscher’s function $`\varphi (k)`$ for this new, less symmetrical case. This has been done without difficulty following the method presented in the original paper, Ref. .

## 3 G-parity boundary condition

It is not immediately obvious how to realize this boundary condition for the pion, because we have control only on the underlying quark fields. However, the G-parity operation gives the solution since under G-parity: $`G|\pi >=|\pi >`$. From this, we can see that anti-periodic boundary conditions for the pion can be achieved by performing a G-parity operation at the boundary—an operation that also can be defined at the quark level: $`u\overline{d}`$ and $`d\overline{u}`$. (Recall that G is defined as the product of charge-conjugation and an isospin rotation, $`G=C\mathrm{exp}i\pi I_y`$.) These unusual boundary were originally proposed and implemented by Wiese and similar boundary conditions have also been used in Ref. .

Since we must also treat the Kaon, we must impose a charge conjugate boundary condition on the strange quark for consistency. This is easily done by including a fictitious charmed quark, degenerate with the strange quark and extending the usual G-parity by treating the $`(c,s)`$ quark pair as an independent iso-doublet. Our initial finite volume eigenstate can then be constructed as a G-even mixture $`(K^0+D^0)/\sqrt{2}`$. However, because our weak decay operator and our final states have definite strangeness and charm, this mixing only introduces an extra “finite-volume factor” of $`1/\sqrt{2}`$ which can easily be removed.

## 4 Lattice Action

### 4.1 Symmetry

By imposing this G-parity boundary condition, we break chiral symmetry in a manner similar to the chiral symmetry breaking of a mass term. This was studied for the continuum in Ref. . On the lattice, the boundary condition is embedded in the action itself, adding a boundary-specific term. For example, with G-parity boundary conditions we have a term in the action density of the form:

$$\overline{u}(L_z1)U^3(L_z1)\gamma _3C\overline{d}(0)^T,$$
(1)
for Dirac conventions in which charge conjugation is represented by $`qC\overline{q}^T`$. Here $`\overline{u}(L_z1)`$ is evaluated at the site $`z=L_z1`$ and $`\overline{d}(0)`$ at the site $`z=0`$ while $`U^3(L_z1)`$ is the link variable joining these two sites across the boundary. Under chiral $`SU(2)_LSU(2)_R`$ such a term transforms as a $`(\frac{1}{2},\frac{1}{2})`$. While such a boundary term may appear to break translational invariance, the fact that the lattice Lagrangian is invariant under the G-parity operation ensures us that the system has translational symmetry, if appropriately defined. (Note, when exploring unconventional boundary conditions it is important to preserve the conservation of lattice momentum and hence to require translational invariance.)

### 4.2 Implementation

The actual implementation of these boundary conditions represents a fairly minor modification of the usual lattice Dirac operator. An easy way of visualizing these G-parity boundary conditions begins with doubling the lattice in the $`z`$-direction, with the complex conjugates of gauge links in the original volume copied onto the links translated by $`L_z`$ into this extended volume. One then easily defines a single-flavor Dirac operator on this $`L_x\times L_y\times 2L_z`$ lattice in the usual way, but imposing anti-periodic boundary conditions between the $`z=0`$ and $`z=2L_z1`$ boundaries. If a quark field in the original volume is identified as a $`u`$ quark, then the quark field in the doubled, $`L_zz<2L_z1`$ volume would be viewed as a $`\overline{d}`$ quark. Note, the computational load is twice that of a usual calculation on a $`L_x\times L_y\times L_z`$ volume—our boundary conditions have effectively required separate treatment for the $`u`$ and $`d`$ quarks. Such an approach works for both Wilson and domain wall fermions.

## 5 Conclusion

Since the G-parity boundary condition is imposed in only one direction, we lose cubic symmetry. In the cubically symmetric case the state of interest with orbital angular momenta $`l=0`$ is mixed only with states having $`l`$ differing by multiples of 4. However, without this cubic symmetry only parity invariance remains which permits mixing with all even $`l`$. Thus, in the practical application of the Lellouch and Lüscher method where the mixed, higher-$`l`$ states are treated as free and their phase shifts ignored, we must neglect the $`l=2`$ phase shift $`\delta _2(m_K)`$. This reduces the accuracy of this approximation, introducing errors of order $`(R/L)^4`$ rather than $`(R/L)^8`$ as is the case for cubic symmetry. Here $`R`$ is a typical strong interaction distance, e.g. $`R=1/m_\rho `$.

There is an additional complexity for the $`I=0`$ final state that should be discussed. For simplicity, we consider a full QCD calculation which includes only dynamical $`u`$ and $`d`$ quarks. For such a system, in addition to the three pions, there will be a fourth $`SU(2)`$ singlet, quark-anti-quark state called here the $`\eta ^{}`$. This $`\eta ^{}`$ state is G-parity even in contrast to G-parity odd triplet of pions and hence will obey even boundary conditions. Thus an intended two-pion, $`I=0`$ state may mix with an unwanted, threshold two-$`\eta ^{}`$ state which will have the same quantum numbers. Fortunately, in the two-flavor calculation just described, axial anomaly effects are expected make the non-Goldstone $`\eta ^{}`$ heavier than the physical Kaon so this two-$`\eta ^{}`$ state will play no role. However, in a quenched calculation, this two-$`\eta ^{}`$ state will be essentially degenerate with the desired $`I=0`$ two-pion state and obstruct the calculation. In fact, the hairpin diagrams entering the quenched two-$`\eta ^{}`$ amplitudes will give an extra power of $`t`$ or $`t^2`$ further obscuring the state of interest. Thus, this G-parity boundary condition can be used in the quenched approximation only for the evaluation of $`\mathrm{\Delta }I=3/2`$ amplitudes.

In summary, we propose imposing anti-periodic boundary conditions on the pion by implementing a G-parity operation at boundary. This reduces the required simulation volume to $`(3\mathrm{Fermi})^3`$ and, more importantly, allows the physical $`\pi \pi (I)|K_{\mathrm{weak}}|K`$ amplitude to be extracted from the leading/next-leading large-time behavior of a Euclidean correlation function for the $`I=2/I=0`$ final states. We hope to try a calculation according to this method in the coming year and thank our colleagues in the RBC collaboration for useful suggestions and criticisms.


format: text

text:
# Anomalous Longitudinal Magnetic Field near the Surface of Copper Conductors

## References


text:
# A New Search Technique for Short Orbital Period Binary Pulsars

## 1 Introduction

Short orbital-period binary pulsars have proven to be sensitive laboratories that allow tests of a variety of physical processes including theories of gravitation. Unfortunately, except for those systems containing the brightest pulsars, they are notoriously difficult to detect. When the observation time $`T_{obs}`$ is greater than a small fraction of the orbital period $`P_{orb}`$, Doppler effects due to the orbital motion of the pulsar cause drastic reductions in the sensitivities of pulsar searches (Johnston & Kulkarni 1991).

Numerous groups (e.g. Middleditch & Kristian 1984; Anderson et al. 1990; Wood et al. 1991) have developed variations of so called “acceleration” searches in order to mitigate the loss in sensitivity caused by orbital motion. These searches can almost completely recover this lost sensitivity when $`T_{obs}P_{orb}/10`$ by taking advantage of the fact that orbital motion causes an approximately linear drift of the pulsar spin frequency $`f_{spin}`$. Recent acceleration searches of globular clusters have discovered numerous binary systems, some with orbital periods as short as $`90`$ min (Camilo et al. 2000; Ransom et al. 2001; D’Amico et al. 2001).

Acceleration searches are effective only for $`T_{obs}P_{orb}/10`$, but search sensitivities improve as $`T_{obs}^{1/2}`$. In effect, acceleration searches are limited to discovering only the brightest binary pulsars in ultra-compact — systems with orbital periods less than a few hours — orbits. For a 2 ms pulsar in a 1 hour orbit with a low mass white dwarf (WD) companion, the optimal observation time $`T_{best}`$ for an acceleration search (using eqn. 22 from Johnston & Kulkarni 1991) is only $`300600`$ seconds. If this pulsar were located in a globular cluster, where observations of $`510`$ hours in duration are common, acceleration searches would be a factor of $`\sqrt{T_{obs}/T_{best}}8`$ times less sensitive than an optimal (i.e. coherent) search of a full observation.

We have developed an efficient search technique for short orbital period binary pulsars that complements the use of acceleration searches. Since it requires $`T_{obs}1.5P_{orb}`$ and its sensitivity improves as the ratio $`r_{orb}=T_{obs}/P_{orb}`$ *increases*, it allows the search of very interesting portions of orbital parameter space for the first time. This technique was developed and tested over the past several years as part of a Ph.D. thesis (Ransom 2001). A similar but independent derivation of many of the properties of this technique has been recently published by Jouteux et al. (2002), based on an earlier mention of this work in Ransom (2000). While the basic search techniques presented by these papers are nearly identical, aspects concerning the determination of the orbital elements of a candidate pulsar discovered with this technique are presented in more detail in this work.

## 2 Phase Modulation due to Orbital Motion

During observations of a binary pulsar where $`T_{obs}`$ is longer than $`P_{orb}`$, the orbital motion causes a modulation of the phase of the pulsar signal. Differences in light travel time across the projected orbit advance or delay the pulse phase in a periodic fashion.

### 2.1 Circular Orbits

If the orbit is circular (where the eccentricity $`e`$ and the angle of periapsis $`\omega `$ are both zero), the phase delays are sinusoidal. The fundamental harmonic111The following analysis can be applied to any harmonic $`h`$ of a pulsar’s signal by substituting $`hr_{spin}`$ for $`r_{spin}`$ and the appropriate harmonic phase for $`\varphi _{spin}`$. of the resulting signal as sampled at a telescope can be described by

$$n_j=n(u)=a\mathrm{cos}\left[2\pi r_{spin}u+\varphi _{spin}+\mathrm{\Phi }_{orb}\mathrm{cos}(2\pi r_{orb}u+\varphi _{orb})\right],$$
(1)
where $`a`$ and $`\varphi _{spin}`$ are the amplitude and phase of the pulsation, $`\mathrm{\Phi }_{orb}`$ and $`\varphi _{orb}`$ are the amplitude and phase of the modulation due to the orbit, $`r`$ is a Fourier frequency (i.e. $`r=fT_{obs}=T_{obs}/P`$), and $`u=j/N=(t_jt_o)/T_{obs}`$ is the discrete dimensionless time during the observation of the $`j^{th}`$ sample of $`N`$ such that $`0u1`$. We can write $`\mathrm{\Phi }_{orb}`$ and $`\varphi _{orb}`$ in terms of the three Keplerian orbital parameters for circular orbits ($`P_{orb}`$, the semi-major axis $`x_{orb}=a_1\mathrm{sin}(i)/c`$, and the time of periapsis $`T_o`$) as

$$\mathrm{\Phi }_{orb}=2\pi x_{orb}f_{spin}7.375\frac{f_1^{1/3}P_{orb}^{2/3}(\mathrm{hr})}{P_{spin}(\mathrm{sec})}$$
(2)
and

$$\varphi _{orb}=2\pi \frac{t_oT_o}{P_{orb}}+\frac{\pi }{2},$$
(3)
where the units of both are radians, and the pulsar’s mass function $`f_1`$ is defined as

$$f_1=\frac{\left[M_2\mathrm{sin}(i)\right]^3}{\left(M_1+M_2\right)^2},$$
(4)
where $`i`$ is the orbital inclination and $`M_1`$ and $`M_2`$ are the pulsar and companion masses (in solar masses) respectfully.

The Discrete Fourier Transform (DFT) of the signal in eqn. 1 when $`\mathrm{\Phi }_{orb}=0`$ can be approximated as (see e.g. Ransom, Eikenberry, & Middleditch 2002)

$`A(r)`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{N1}{}}}n_j\mathrm{e}^{2\pi \mathrm{i}jr/N}N{\displaystyle _0^1}n(u)\mathrm{e}^{2\pi \mathrm{i}ru}du`$ (5a)
$``$ $`{\displaystyle \frac{aN}{2}}\mathrm{e}^{\mathrm{i}\varphi _{spin}}{\displaystyle _0^1}\mathrm{e}^{2\pi \mathrm{i}(r_{spin}r)u}du`$ (5b)
$`=`$ $`A_{coherent}\mathrm{e}^{\mathrm{i}\pi (r_{spin}r)}{\displaystyle \frac{\mathrm{sin}\left[\pi (r_{spin}r)\right]}{\pi (r_{spin}r)}},`$ (5c)
where $`A_{coherent}=\frac{aN}{2}\mathrm{e}^{\mathrm{i}\varphi _{spin}}`$ represents the coherent Fourier response at the spin frequency $`r_{spin}`$. Similarly, we can approximate the DFT of our phase modulated signal from eqn. 1 as

$`A(r)`$ $`=`$ $`{\displaystyle \frac{aN}{2}}{\displaystyle _0^1}\mathrm{e}^{\mathrm{i}\left[2\pi r_{spin}u+\varphi _{spin}+\mathrm{\Phi }_{orb}\mathrm{cos}(2\pi r_{orb}u+\varphi _{orb})\right]}\mathrm{e}^{2\pi \mathrm{i}ru}du`$ (6a)
$`=`$ $`A_{coherent}{\displaystyle _0^1}\mathrm{e}^{\mathrm{i}\mathrm{\Phi }_{orb}\mathrm{cos}(2\pi r_{orb}u+\varphi _{orb})}\mathrm{e}^{2\pi \mathrm{i}(r_{spin}r)u}du.`$ (6b)
Using the Jacobi-Anger expansion (e.g. Arfken 1985)

$$\mathrm{e}^{\mathrm{i}z\mathrm{cos}(\theta )}=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{i}^mJ_m(z)\mathrm{e}^{\mathrm{i}m\theta }=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}J_m(z)\mathrm{e}^{\mathrm{i}m\left(\theta +\frac{\pi }{2}\right)},$$
(7)
where $`J_m(z)`$ is an integer order Bessel function of the first kind, we can expand the integrand in eqn. 6b to give

$`A(r)`$ $`=`$ $`A_{coherent}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}J_m(\mathrm{\Phi }_{orb}){\displaystyle _0^1}\mathrm{e}^{\mathrm{i}m(2\pi r_{orb}u+\varphi _{orb}+\frac{\pi }{2})}\mathrm{e}^{2\pi \mathrm{i}(r_{spin}r)u}du`$ (8a)
$`=`$ $`A_{coherent}{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}J_m(\mathrm{\Phi }_{orb})\mathrm{e}^{\mathrm{i}m\left(\varphi _{orb}+\frac{\pi }{2}\right)}{\displaystyle _0^1}\mathrm{e}^{2\pi \mathrm{i}\left(r_{spin}+mr_{orb}r\right)u}du.`$ (8b)
Since the integrand in eqn. 8b is identical in form to that in eqn. 5b, the Fourier transform of a phase modulated signal is equivalent to the Fourier transform of a series of cosinusoids at frequencies centered on the pulsar spin frequency $`r_{spin}`$, but separated from $`r_{spin}`$ by $`mr_{orb}`$ Fourier bins (i.e. sidebands). The response is therefore

$$A(r)=A_{coherent}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}J_m(\mathrm{\Phi }_{orb})\mathrm{e}^{\mathrm{i}m\left(\varphi _{orb}+\frac{\pi }{2}\right)}\mathrm{e}^{\mathrm{i}\pi (r_{spin}+mr_{orb}r)}\frac{\mathrm{sin}\left[\pi (r_{spin}+mr_{orb}r)\right]}{\pi (r_{spin}+mr_{orb}r)}.$$
(9)
When $`r_{orb}`$ is an integer (i.e. the observation covers an integer number of complete orbits), the Fourier response shown in eqn. 9 becomes particularly simple at Fourier frequencies $`r_{spin}+sr_{orb}`$, where $`s`$ is an integer describing the sideband in question. In this case, the summation collapses to a single term when $`m=s`$, giving

$$A\left(r_{spin}+sr_{orb}\right)=A_{coherent}J_s(\mathrm{\Phi }_{orb})\mathrm{e}^{\mathrm{i}s\left(\varphi _{orb}+\frac{\pi }{2}\right)}.$$
(10)
In effect, the phase modulated signal produces sidebands composed of a series of sidebands split in frequency from $`r_{spin}`$ by $`sr_{orb}`$ Fourier bins, with amplitudes proportional to $`J_s(\mathrm{\Phi }_{orb})`$, and phases of $`s\left(\varphi _{orb}+\frac{\pi }{2}\right)+\varphi _{spin}`$ radians222Note that $`J_n(z)=(1)^nJ_n(z)`$, and also that $`J_n(z)`$ can be negative. If the sideband amplitude determined by the Bessel function is negative (i.e. $`J_s(z)`$ is negative), the *measured* phase of the Fourier response will differ by $`\pi `$ radians from the predicted phase since the *measured* amplitude is always positive.. Figure 1 shows the Fourier response for a “typical” pulsar-WD system in a circular 50 min orbit during an 8 hour observation.

While the number of sidebands implied by eqn. 9 is infinite (assuming $`\mathrm{\Phi }_{orb}0`$), properties of the Bessel functions — which define the magnitude of the Fourier amplitudes of the sidebands — give a rather sharp cutoff at $`2\mathrm{\Phi }_{orb}`$ sidebands (or $`\mathrm{\Phi }_{orb}`$ pairs of sidebands). The maximum value of $`J_s(\mathrm{\Phi }_{orb})`$ also occurs near $`s=\mathrm{\Phi }_{orb}`$ at a value of $`0.675\mathrm{\Phi }_{orb}^{1/3}`$, which gives the phase modulation response its distinctive “horned” shape when $`\mathrm{\Phi }_{orb}`$ is large (see Fig. 1). The average magnitude of the phase modulation sidebands falls off more quickly as $`0.611\mathrm{\Phi }_{orb}^{1/2}`$ (see Fig. 2).

For pulsars with narrow pulse shapes and many Fourier harmonics, phase modulation produces sidebands split by $`r_{orb}`$ Fourier bins centered around each spin harmonic. For spin harmonic number $`h`$, the phase modulation amplitude is $`\mathrm{\Phi }_{orb,h}=2\pi x_{orb}hf_{spin}`$. Therefore, for large $`\mathrm{\Phi }_{orb,h}`$, there are $`h`$ times more sidebands than around the fundamental, but the average amplitudes are smaller by a factor of $`h^{1/2}`$.

### 2.2 Small Amplitude Limit

When $`\mathrm{\Phi }_{orb}1`$, a phase modulated signal will have only three significant peaks in its Fourier response. These peaks are located at $`r_{spin}`$ and $`r_{spin}\pm r_{orb}`$ and have magnitudes proportional to $`J_0(\mathrm{\Phi }_{orb})1\mathrm{\Phi }_{orb}^2/41`$ and $`J_{\pm 1}(\mathrm{\Phi }_{orb})\pm \mathrm{\Phi }_{orb}/2`$ respectively.

Middleditch et al. (1981) performed a detailed analysis of phase modulated optical pulsations in the small amplitude limit from the 7.7 second X-ray pulsar 4U 1626-67 in order to solve for the orbital parameters of the system. The separation of the sidebands provided a measurement of the $`42`$ min orbital period and the magnitudes and phases of the three significant Fourier peaks allowed them to solve for $`\mathrm{\Phi }_{orb}`$ (and therefore $`x_{orb}`$) and $`\varphi _{orb}`$ directly.

Small amplitude phase modulation might also be observed from freely precessing neutron stars (e.g. Nelson, Finn, & Wasserman 1990; Glendenning 1995). For certain geometries, the geometric wobble of the pulsar beam would periodically modulate the pulse arrival times and cause the system to appear as a compact binary. If such rapidly precessing objects exist, they would display a small number ($`13`$) of pairs of sidebands in a power spectrum of a long observation, since the maximum possible phase modulation amplitude for such a system is $`\mathrm{\Phi }_{orb}\pi /2`$ radians (Nelson et al. 1990).

### 2.3 Elliptical Orbits

Konacki & Maciejewski (1996) derived a Fourier expansion for elliptical orbital motion with the useful property that the magnitudes of higher order harmonics of the orbital expansion (which correspond to independent modulation amplitudes $`\mathrm{\Phi }_{orb,h^{}}`$) decrease monotonically. In the limit of circular orbits the expansion reduces to the single phase modulating sinusoid found in eqn. 1.

Each harmonic of the orbital modulation generates sidebands around the spin harmonic of the form discussed in §2.1. The sidebands created by harmonic $`h^{}`$ of the orbital expansion are separated from each other by $`h^{}r_{orb}`$ Fourier bins. Since $`\mathrm{\Phi }_{orb,h^{}}`$ monotonically decreases with increasing orbital harmonic $`h^{}`$, the number of sidebands generated decreases while their spacing and average amplitudes increase. The most important aspect of this superposition of sidebands is the fact that the Fourier response near a spin harmonic contains sidebands spaced by $`r_{orb}`$ Fourier bins — from the fundamental of the orbital expansion — just as in the limiting case of circular orbits.

## 3 Pulsars with Deterministic Periodic Phase Modulation

Figures 3 and 4 show $`P_{spin}`$, $`P_{orb}`$, and $`\mathrm{\Phi }_{orb}`$ for 54 binary pulsars from the literature, all with orbital periods less than 10 days. A few conclusions can be reached immediately. First, radio millisecond pulsars (MSPs) are the most common binary pulsar. Second, there is a relatively flat distribution of $`P_{orb}`$ for radio pulsars which cuts off quite dramatically near $`P_{orb}0.08`$ d ($`2`$ h) — the two systems with shorter orbital periods are the recently discovered accreting X-ray MSPs XTE J1751-314 (Markwardt et al. 2002) and XTE J0929-314 (Galloway et al. 2002). And third, a large fraction of the systems are located in globular clusters.

The fact that large numbers of relatively compact binary pulsar systems should contain MSPs is not surprising in light of the standard “recycling” model for MSP creation (Verbunt 1993; Phinney & Kulkarni 1994). This model explains not only the preponderance of pulsars with millisecond spin periods, but also the fact that many MSPs have low-mass ($`0.020.4\text{M}_{}`$) companion stars in compact circular orbits. These recycled systems make up the vast majority of the pulsars shown in Figure 3 and have typical values of $`\mathrm{\Phi }_{orb}`$ of $`10^2`$ to a few times $`10^3`$ radians.

The cutoff in $`P_{orb}`$ at approximately two hours is almost certainly due in some part to selection effects. The rapid orbital motion of these systems prevents their detection using conventional search techniques — including acceleration searches — for all but the brightest systems (Johnston & Kulkarni 1991). In fact, the only MSPs discovered with shorter orbital periods ($`P_{orb}40`$ min for both XTE J1751-314 and XTE J0929-314) were discovered during transient X-ray outbursts when the pulsed fractions were very high (Markwardt et al. 2002; Galloway et al. 2002). Since selection effects keep us from detecting ultra-compact systems that are not exceptionally bright, it is difficult to determine how much the observed lack of such systems is due to a true lack of compact systems. It is certain that neutron star (NS) systems with much shorter orbital periods *do* exist, as the new accreting X-ray MSPs and other ultra-compact LMXBs such as 4U 1820$``$303 ($`P_{orb}11`$ min), 4U 1850$``$087 ($`P_{orb}21`$ min), 4U 1627$``$673 ($`P_{orb}42`$ min), and 4U 1916$``$053 ($`P_{orb}50`$ min) (Chou et al. 2001; Liu et al. 2001) have shown.

The large fraction of known systems inhabiting globular clusters can be explained by three basic facts. First, due to their known location on the sky, long targeted observations of clusters allow fainter systems to be detected (e.g. Anderson 1992). Second, repeated observations of certain clusters have allowed many very weak pulsars to be detected with the help of amplification from interstellar scintillation (Camilo et al. 2000; Possenti et al. 2001). Finally, large numbers of binary pulsars are expected to have been created in globular clusters due to dynamical interactions (e.g. Kulkarni & Anderson 1996; Davies & Hansen 1998; Rasio et al. 2000).

Recently, Rasio et al. (2000) have conducted extensive population synthesis studies for binary pulsars in dense globular clusters like 47 Tucanae. Their initial results imply that exchange interactions among primordial binaries will produce large numbers of NSs in binary systems. In fact, their results predict that large populations of binary pulsars with low mass companions ($`m_c0.05\text{M}_{}`$) should exist in many globular clusters with orbital periods as short as $`15`$ min.

The evidence seems to suggest that a large population of ultra-compact binary pulsars should exist — particularly in globular clusters — which have so far eluded detection. With long observations, these pulsars will be detected if search algorithms can identify the distinctive sidebands generated by the systems. The rest of this paper discusses such an algorithm.

## 4 Sideband Searches

The fact that power from an orbitally modulated sinusoid is split into multiple sidebands allows us to increase the signal-to-noise of the detected response by incoherently summing the sidebands. According to one of the addition theorems of Bessel functions

$$\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}J_m(z)^2=1,$$
(11)
theoretically *all* of the signal power can be recovered in this manner for noiseless data. Practically, though, complete recovery of signal power is not possible for data with a finite signal-to-noise ratio — especially if $`\mathrm{\Phi }_{orb}1`$. Such a summation would require calculation of the $`2\mathrm{\Phi }_{orb}`$ significant sidebands at the precise Fourier frequencies of their peaks $`r_{spin}\pm sr_{orb}`$. This computationally daunting task would result in very large numbers of independent search trials when $`\mathrm{\Phi }_{orb}`$, $`r_{spin}`$, and $`r_{orb}`$ are unknown.

Even if we could recover all of the power using a sum of the sidebands, the significance of the measurement would be less than the significance of a non-modulated sinusoid with the same Fourier power. This loss in significance is due to the fact that noise is co-added along with signal. The exact loss in significance can easily be calculated since the probability for a sum of $`m`$ noise powers to exceed some power $`P_{sum}`$ is the probability for a $`\chi ^2`$ distribution of $`2m`$ degrees of freedom to exceed $`2P_{sum}`$ (see e.g. Groth 1975, and references therein).

In the case where only a few sidebands are suspected to be present (i.e. $`\mathrm{\Phi }_{orb}5`$), a brute-force search using incoherent sideband summing may be computationally feasible and worthwhile. Unfortunately, except for binary X-ray pulsars with long spin periods and possibly certain freely precessing NSs (see §2.2), the number of astrophysical systems with $`\mathrm{\Phi }_{orb}5`$ is probably small compared to those with $`\mathrm{\Phi }_{orb}1`$ (see §3).

### 4.1 Two-Staged Fourier Analysis

For systems with $`\mathrm{\Phi }_{orb}1`$, the large number of regularly spaced sidebands allows a very efficient detection scheme. *Since the sidebands are separated from each other by $`r_{orb}`$ Fourier bins, they appear as a localized periodicity in the power spectrum of the original time series which can be detected with a second stage of Fourier analysis.* By stepping through the full length power spectrum and taking short Fast Fourier Transforms (FFTs) of the powers, we incoherently sum any sidebands present and efficiently detect their periodic nature. As an added benefit, the measured “frequency” of these “sideband pulsations” gives a direct and accurate measurement of the *orbital* period.

We define the following notation to represent the two distinct stages of the Fourier analysis for detecting a phase modulated signal. As defined in eqn. 1, the initial time series $`n_j`$ contains $`N`$ points. After Fourier transforming the $`n_j`$ (usually using an FFT of length $`N`$), the complex responses at Fourier frequencies $`r`$ are represented by $`A(r)`$ as defined in eqn. 5a. The powers and phases are then simply $`P(r)`$ and $`\varphi (r)`$ respectively. The short FFTs of the $`P(r)`$ used to detect the sideband pulsations (i.e. the secondary power spectrum) are of length $`M`$ (where $`MN`$ and should closely match the extent of the sideband structure in the original power spectrum) and generate complex responses $`A_2(r)`$, powers $`P_2(r)`$, and phases $`\varphi _2(r)`$.

In order for a pulsed signal to undergo enough modulation to produce periodic sidebands, the observation must contain more than one complete orbit (i.e. $`r_{orb}>1`$). One might initially think that two complete orbits would be necessary to create a Nyquist sampled series of sidebands in the original power spectrum. However, since the sideband pulsations are not bandwidth limited, periodic signals with “wavelengths” less than two Fourier bins (i.e. $`r_{orb}<2`$) will still appear in the $`P_2(r)`$. Instead of appearing at Fourier frequency $`r=M/r_{orb}`$ as for signals with $`r_{orb}>2`$ (see §4.2.2), the fundamental harmonic of the response will appear aliased around the Nyquist frequency ($`r_{Nyq}=M/2`$) at $`r_{alias}=MM/r_{orb}`$. While $`r_{orb}>1`$ is required to begin to produce periodic sidebands, simulations have shown that this technique realistically requires $`r_{orb}1.5`$ to create detectable sidebands from most binary systems during reasonable observations.

As the number of orbits present in the data ($`r_{orb}`$) increases, so does the spacing (in Fourier bins) between the sidebands. Since each sideband has a traditional sinc function shape in $`A_2(r)`$ (see §2.1), the full-width at half-maximum (FWHM) of a sideband is approximately one Fourier bin. Therefore, an increase in spacing effectively decreases the duty-cycle of the sideband pulsations (i.e. $`1/r_{orb}`$). Signals with small duty-cycles are easier to detect for two reasons. First, significant higher harmonics begin to appear in $`P_2(r)`$ which can be incoherently summed in order to improve the signal-to-noise ratio. A rule-of-thumb for the number of higher harmonics present in $`P_2(r)`$ is approximately one-half the inverse of the duty-cycle ($`r_{orb}/2`$). Second, as the duty-cycle decreases, the Fourier amplitudes of each harmonic increase in magnitude up to a maximum of twice that of a sinusoid with an identical pulsed fraction (see Ransom et al. 2002, for a more detailed discussion of pulse duty-cycle effects). This effect makes each individual harmonic easier to detect as $`r_{orb}`$ increases.

### 4.2 Search Considerations

As the number of orbits present in a time series increases, the signal-to-noise ratio of a detectable pulsar spin harmonic goes approximately as $`(r_{orb}1)^\beta `$ where $`0.3\beta 0.4`$ (see Figures 6 and 7). In contrast, the signal-to-noise ratio of a pulsar spin harmonic detected coherently goes as $`T_{obs}^{1/2}`$. This leads to a rule-of-thumb for those conducting sideband searches: “The longer the observation the better.” When making very long observations it is important to keep the on-source time (i.e. the window function) as continuous and uniform as possible. Significant gaps in the data — especially periodic or recurring gaps of the kind common in X-ray observations — create side-lobes around Fourier peaks of constant-frequency pulsations that a sideband search will detect. It is a good idea before a sideband search to take the power spectrum of the power spectrum of the known window function in order to identify spurious orbital periods that a sideband search might uncover.

Since phase modulation sidebands require coherent pulsations (i.e. $`f_{spin}`$ is constant in time) in order to form, it may be necessary to correct the time series (or equivalently the Fourier transform) such that the underlying pulsations are coherent. As an example, if the telescope’s motion with respect to the solar system barycenter is not taken into account during long observations, interstellar pulsations — and therefore any phase modulation sidebands as well — will be “smeared” across numerous Fourier bins due to Doppler effects. By stretching or compressing the time series to account for the Earth’s motion before Fourier transforming (or equivalently, by applying the appropriate Fourier domain matched filter *after* the Fourier transform, see Ransom et al. 2002) the smeared signal can be made coherent for the purposes of the sideband search.

If the time-dependence of $`f_{spin}`$ is unknown *a priori*, due to accretion in an LMXB, the spin-down of a young pulsar, or timing (spin) noise, for example, one could in principle attempt a series of trial corrections and a sideband search for each trial. Such a search methodology is similar in computational complexity to traditional “acceleration” searches.

While §4.1 provided basic principles for detecting compact binary pulsars using the two-staged Fourier analysis, in the rest of this section we provide a more detailed discussion of how a search might be conducted. Figure 5 shows two examples of sideband searches. The left column shows data from an 8 h observation of the globular cluster 47 Tucanae taken with the Parkes radio telescope at 20 cm on 2000 November 17. The signal is the fundamental harmonic of the 2.1 ms binary pulsar 47 Tuc J ($`f_{spin}=476.05`$ Hz, $`P_{orb}=2.896`$ h, $`r_{orb}=2.762`$, $`\mathrm{\Phi }_{orb}=120.87`$ radians, and $`m_c0.03`$$`\text{M}_{}`$) with a signal-to-noise ratio of $`\eta =a/\sigma _{n_j}0.0054`$. The right column shows data from a simulated 8 h observation similar to the 47 Tuc observation described above. The signal is the fundamental harmonic of the 2.0 ms binary pulsar used in Figure 1 ($`f_{spin}=500`$ Hz, $`P_{orb}=50`$ min, $`r_{orb}=8`$, $`\mathrm{\Phi }_{orb}=413.44`$ radians, and $`m_c=0.2`$$`\text{M}_{}`$) with a signal-to-noise ratio one half that of the 47 Tuc J observation ($`\eta 0.0027`$).

Once a long observation has been obtained and prepared as discussed above, we Fourier transform the data and compute the power spectrum $`P(r)`$. We then begin a sideband search which is composed of two distinct parts: initial detection of a binary pulsar and the determination of a detected pulsar’s orbital elements.

#### 4.2.1 Initial Detection

To maximize the signal-to-noise of a detection in $`P_2(r)`$, we want to match both the length and location of a short FFT $`M`$ with the width and location of a pulsar’s sideband response ($`2r_{orb}\mathrm{\Phi }_{orb}`$) in $`P(r)`$. Because a pulsar’s orbital parameters are unknown *a priori*, we must search a range of short FFT lengths $`M`$ and overlap consecutive FFTs. Therefore, we must choose the number and lengths of the short FFTs, the fraction of each FFT to overlap with the next, and the number of harmonics to sum in each FFT. These choices are influenced by the length of the observation and the nature of the systems that we are attempting to find.

All currently known binary pulsars with $`P_{orb}<10`$ d have modulation amplitudes of $`\mathrm{\Phi }_{orb}6500`$ radians (see Figure 4), with the vast majority in the range $`\mathrm{\Phi }_{orb}205000`$ radians. This range also includes many currently unknown but predicted “holy grail”-type systems such as compact MSP-black hole systems or MSP-WD systems with orbital periods of $`10`$ min. This range of $`\mathrm{\Phi }_{orb}`$ values is due to the fact that the short spin periods of MSPs tend to offset the small semi-major axes ($`x_{orb}`$) of their typical binary systems, while binary systems with wider orbits tend to have longer spin periods. Multiplying this range by twice the minimum $`r_{orb}`$ of 1.5 gives a rule-of-thumb range of $`M6015000`$ where the high end should be increased based on expected values of $`r_{orb}`$. For $`812`$ h globular cluster searches a reasonable range of values would be $`64M65536`$ in powers-of-two increments, for a total of 11 different values for $`M`$.

The choice of the overlap fraction is more subjective. To optimize a detection in $`P_2(r)`$, one should overlap the short FFTs by at least 50% in order to avoid excessive loss of signal-to-noise when a series of sidebands straddles consecutive FFTs. Larger overlap percentages may result in higher signal-to-noise, but at the expense of increased computational costs and numbers of search trials (even though the trials are not completely independent when overlapping).

In Figure 5 plots $`a`$ and $`e`$, the grey shaded regions show the sections of the $`P(r)`$ that were Fourier transformed for the sideband search. The short FFTs of length $`M=1024`$ and $`M=8192`$ were centered on the known pulsar spin frequency and resulted in the $`P_2(r)`$ displayed in plots $`b`$ and $`f`$ respectively. The short FFT lengths were chosen as the shortest powers-of-two greater than $`2r_{orb}\mathrm{\Phi }_{orb}`$. For 47 Tuc J (plot $`a`$), the “horned” structure of the phase modulation sidebands described in §2.1 is easily recognizable. This shape is not apparent in plot $`e`$, and in fact, traditional searches would detect nothing unusual about this section of the power spectrum.

Normalization of the $`P_2(r)`$ results in the usual exponential distribution with mean and standard deviation of one for a transform of pure noise. If one first normalizes the $`P(r)`$ as shown in plots $`a`$ and $`e`$ (see Ransom et al. 2002, for several normalization techniques), the resulting $`P_2(r)`$ can be normalized simply by dividing by $`M`$.

Plots $`b`$ and $`f`$ of Figure 5 show the highly significant detections of the binary pulsars in $`P_2(r)`$. For 47 Tuc J, only one orbital harmonic is visible in $`P_2(r)`$ since $`r_{orb}<2`$, but the single-trial significance is approximately $`23\sigma `$. For the simulated 50 min binary pulsar, 3 significant orbital harmonics are detected (predicted harmonic locations are marked with a dotted line in plot $`f`$). When summed, the detection has a single-trial significance of $`15\sigma `$.

#### 4.2.2 Orbit Determination

Once a binary pulsar has been detected, examination of the $`P(r)`$ and $`\varphi (r)`$ allow the estimation of the Keplerian orbital parameters and the pulsar spin period $`P_{spin}`$. Brute-force matched filtering searches centered on the estimates are then used to refine the parameters and recover the fully coherent response of the pulsar.

The simplest parameter to determine is the orbital period $`P_{orb}`$. The spacing between sidebands is $`r_{orb}=T_{obs}/P_{orb}`$ bins, implying that the most significant peak in $`P_2(r)`$ should occur at Fourier frequency $`M/r_{orb}`$. Conversely, if we measure the Fourier frequency $`r_{meas}`$ in $`P_2(r)`$ with the most power, we can solve for $`r_{orb}`$ using $`r_{orb}=M/r_{meas}`$. By using Fourier interpolation or zero-padding to oversample $`A_2(r)`$ and $`P_2(r)`$, we can determine $`r_{meas}`$ to an accuracy of

$$\sigma _r=\frac{3}{\pi \alpha \sqrt{6P_2(r_{meas})}},$$
(12)
where $`P_2(r_{meas})`$ is the normalized power at the peak and $`\alpha `$ is the signal “purity” — a property that is proportional to the root-mean-squared (RMS) dispersion of the pulsations in time about the centroid — which is equal to one for pulsations present throughout the data (see Middleditch et al. 1993; Ransom et al. 2002). Remembering that $`P_{orb}=T_{obs}/r_{orb}`$ and solving gives

$$P_{orb}=\frac{T_{obs}r_{meas}}{M}\pm \frac{3T_{obs}}{\pi \alpha M\sqrt{6P_2(r_{meas})}}.$$
(13)
For aliased signals, the orbital period is

$$P_{orb,aliased}=\frac{T_{obs}\left(Mr_{meas}\right)}{M}\pm \frac{3T_{obs}}{\pi \alpha M\sqrt{6P_2(r_{meas})}}.$$
(14)
For pulsars with $`\mathrm{\Phi }_{orb}1`$ observed such that $`r_{orb}>2`$, it is often possible to measure $`P_{orb}`$ to one part in $`10^4`$, corresponding to a $`1`$ s error for many of the pulsars shown in Figure 4.

Once $`P_{orb}`$ is known it becomes possible to measure the projected semi-major axis $`x_{orb}`$ by determining the total width of the sidebands in $`P(r)`$. For bright pulsars like 47 Tuc J, the number of Fourier bins between the “horns” ($`\mathrm{\Delta }_r2\mathrm{\Phi }_{orb}r_{orb}`$) as shown in Figure 5$`a`$ can be measured directly and converted into $`x_{orb}`$ using

$$x_{orb}\frac{\mathrm{\Delta }_r}{4\pi f_{spin}r_{orb}},$$
(15)
where $`f_{spin}`$ can be estimated as the frequency midway between the phase modulation “horns”.

For weaker pulsars, where the sideband edges are not obvious (as in Figure 5$`e`$), $`\mathrm{\Delta }_r`$ can be estimated using numerous short FFTs. By taking a series of short FFTs of various lengths around the detection region of $`P(r)`$ and measuring the power, “centroid” and “purity”333The measured values of “centroid” and “purity” are estimates of a signal’s location and duration in a time series as determined from the derivatives of the Fourier phase and power at the peak of the signal’s response (see Middleditch et al. 1993; Ransom et al. 2002). at the frequency corresponding to $`P_{orb}`$, one can map the extent of the hidden sidebands and compute an estimate for $`x_{orb}`$.

For pulsars in circular orbits, which constitute the majority of the systems in Figures 3 and 4, the only remaining orbital parameter is the time of periapsis passage $`T_o`$. Defining the time since periapsis passage as $`\mathrm{\Delta }_{T_o}=t_oT_o`$ and using eqn. 3 we get

$$\mathrm{\Delta }_{T_o}=\frac{P_{orb}}{2\pi }\left(\varphi _{orb}\frac{\pi }{2}\right),$$
(16)
where $`\varphi _{orb}`$ can be measured from the phases of the sidebands in $`A(r)`$ using the following technique.

From eqn. 10 we know that the sidebands have phases of $`\varphi _s=s\left(\varphi _{orb}+\frac{\pi }{2}\right)+\varphi _{spin}`$ radians. Unfortunately, some of the measured phases will differ by $`\pi `$ radians from those predicted by eqn. 10 since the corresponding sideband amplitudes (the $`J_s(\mathrm{\Phi }_{orb})`$) are negative (see footnote 5). If we predict values for $`J_s(\mathrm{\Phi }_{orb})`$ using our knowledge of $`r_{orb}`$, $`f_{spin}`$, and $`x_{orb}`$, we can “flip” (i.e. add or subtract $`\pi `$ radians to) the measured phases of the sidebands that are predicted to have negative amplitudes. We can then estimate $`\varphi _{orb}`$ simply by subtracting the phases of neighboring sidebands using $`\varphi _{orb}=\varphi _{s+1}\varphi _s\pi /2`$.

For weak pulsars each measurement of $`\varphi _{orb}`$ has a large uncertainty. Fortunately, we can make $`2\mathrm{\Phi }_{orb}`$ measurements of $`\varphi _{orb}`$ by using each pair of sidebands and then determine $`\varphi _{orb}`$ and its uncertainty statistically. Figure 5 plots $`c`$ and $`g`$ show histograms of the measurement of $`\varphi _{orb}`$ using this technique.

It is important to mention that for weak pulsars it may be difficult to determine the location of the sidebands in order to measure their phases. The location of the sideband peaks can be calculated by measuring the phase $`\varphi _{meas}`$ of the fundamental harmonic discovered in $`A_2(r)`$. The number of Fourier bins from the first bin used for the short FFT to the peak of the first sideband in $`A(r)`$ is simply $`\varphi _{meas}r_{orb}/2\pi `$. Subsequent peaks are located at intervals of $`r_{orb}`$ Fourier bins.

Once estimates have been made for the Keplerian orbital parameters, we compute a set of complex sideband templates (i.e. matched filters for the Fourier domain response of an orbitally modulated sinusoid) over the most likely ranges of the parameter values given the uncertainties in each. These templates are then correlated with a small region of the Fourier amplitudes $`A(r)`$ around $`f_{spin}`$. When a template matches the pulsar sidebands buried in $`A(r)`$, we recover essentially all the power in that pulsar spin harmonic (Ransom et al. 2002). Figure 5 plots $`d`$ and $`h`$ show the results of just such a matched filtering operation. The Doppler effects from the orbital motion have been completely removed from the data and the resulting Fourier response is that of an isolated pulsar.

For binary pulsars in eccentric orbits, the techniques discussed here could in principle be applied to each set of sidebands from the orbital Fourier expansion (see §2.3). Such an analysis would be much more complicated than that described here and would require a high signal-to-noise detection.

## 5 Discussion

We have described a new search technique for binary pulsars that can identify sidebands created by orbital modulation of a pulsar signal when $`T_{obs}>P_{orb}`$. Sideband or phase-modulation searches allow the detection of short orbital period binary pulsars that would be undetectable using conventional search techniques.

Figures 6 and 7 show how the sensitivity of sideband searches compares to that of acceleration searches for a wide range of observation times $`T_{obs}`$. Acceleration search sensitivities are near-optimal when $`T_{obs}0.1P_{orb}`$ (Johnston & Kulkarni 1991), but degrade rapidly at longer integration times when the constant frequency derivative approximation breaks down. Sideband searches, on the other hand, require $`T_{obs}1.5P_{orb}`$ in order to approach the sensitivity of optimal duration acceleration searches, but for longer $`T_{obs}`$ sensitivity improves as $`(r_{orb}1)^\beta `$ where $`0.3\beta 0.4`$. *For targeted searches of duration $``$8$``$12 h (e.g. typical globular cluster observations), sideband searches are $``$2$``$10 times more sensitive to compact binary pulsars ($`P_{orb}4`$ h) than optimal duration acceleration searches.* In general, one should think of acceleration and sideband searches as complementary to each other — sideband searches allow the detection of ultra-compact binary pulsars, while acceleration searches maximize the detectability of isolated and longer-period binary pulsars.

The fact that sideband searches target a different portion of orbital parameter space than acceleration searches and yet require significantly less computation time, provides good reason to include them in future pulsar searches where $`T_{obs}30`$ min. In fact, the Parkes Multibeam Pulsar Survey ($`T_{obs}35`$ min, Lyne et al. 2000; Manchester et al. 2001) has included a sideband search in an on-going re-analysis of their survey data at the cost of only a marginal increase in computer time. While the probability of detecting a binary pulsar with $`P_{orb}23`$ min is almost certainly quite low, the discovery of a single such system would provide a wealth of scientific opportunity.

*Acknowledgments* We would like to thank G. Fazio for supporting our research and encouraging us in this work. Many of the computations for this paper were performed on equipment purchased with NSF grant PHY 9507695. SMR is a Tomlinson Fellow, JMC is supported by NSF Grant 9819931 and by the National Astronomy and Ionosphere Center, and SSE is supported by a NSF CAREER Grant.


text:
# Foam: A General purpose Monte Carlo Cellular AlgorithmWork partly supported by the European Community’s Human Potential Programme under contract HPRN-CT-2000-00149 “Physics at Colliders”, by Polish Government grants KBN 5P03B09320 and 2P03B00122, and by NATO grant PST.CLG.977751.

## Abstract

A general-purpose, self-adapting Monte Carlo (MC) algorithm implemented in the program Foam is described. The high efficiency of the MC, that is small maximum weight or variance of the MC weight is achieved by means of dividing the integration domain into small cells. The cells can be $`n`$-dimensional simplices, hyperrectangles or a Cartesian product of them. The grid of cells, “foam”, is produced in the process of the binary split of the cells. The next cell to be divided and the position/direction of the division hyperplane is chosen by the algorithm which optimizes the ratio of the maximum weight to the average weight or (optionally) the total variance. The algorithm is able to deal, in principle, with an arbitrary pattern of the singularities in the distribution.

Introduction and motivation
For the problem of function minimalization of the arbitrary user-function one applies MINUIT or other program from the NAG library. One finds also many general purpose programs for integration of “arbitrary” integrand. The above general-purpose tools work, in principle, for a very wide range of the user-functions. For multi-dimensional Monte Carlo simulation problem, that is for the problem of generating randomly points according to $`n`$-dimensional distribution, there is precious little examples of the general purpose Monte Carlo Simulators (GPMCS), that is programs which work (in principle) for an arbitrary integrand.

Inevitably the GPMCS has to work in 2 stages: “exploration” and “generation”. During exploration GPMCS is “digesting” the entire shape of the $`n`$-dimensional distribution $`\rho (x_1,x_2,\mathrm{}x_n)`$ to be generated and memorizes it as efficiently as possible using all CPU power and memory available. Obviously, for the memorized $`\rho ^{}(x_1,x_2,\mathrm{}x_n)`$ a method of the MC generation of the points $`\stackrel{}{x}`$ exactly according to $`\rho ^{}(\stackrel{}{x})`$, has to be available. The quality of the distribution of the weight $`w=\rho /\rho ^{}`$ for events in the generation (small variance, good ration of maximum to average, etc.) is determined by the algorithm of the exploration. In other words, the quality of the “target weight distribution” in the latter generation is determining/driving the algorithm of the former exploration.

At the time of writing up this contribution a detailed description of the algorithm and the program is available in ref. . The original, less advanced version was described in ref. . The other old and new efforts in the area of general purpose MC methods see refs.

Cellular algorithm of the Foam
The most obvious method to minimize the variance (or maximum weight) of the target weight distribution in generation is to split integration domain into many cells111This general idea is at least 40 years old., such that $`\rho (\stackrel{}{x})`$ is approximated by constant $`\rho ^{}(\stackrel{}{x})`$ within each cell. This is a cellular class of general purpose MC algorithms222“Stratified sampling”, used in the literature, has a narrower meaning. which always involves a kind of cellular exploration of the user defined distribution.

The obvious questions are: what kind (shape) of cells to use and how to cover the space with cells? In Foam algorithm/program three types of cells are used: pure simplices, pure hyperrectangles and Cartesian products of them. Cells of this kind can be parametrized rather easily and the related computer memory consumption is not excessive.

The system of cells can be created all at once (like in VEGAS ) or in an evolutionary way, by the “splitting process”. The Foam algorithm does the binary split of cells. Next cell to be split is chosen by examining the “target weight distribution”, that is of the event generation in the latter phase. The binary split provides automatically the full coverage of the space, simply because the primary “root cell” is identical with the entire integration domain and there is never any mismatch between parent cell and two daughter cells.

Another important issue is the question of the variance reduction versus maximum weight reduction. In construction of the Foam algorithm we have put most effort on the minimization of the ratio of the maximum weight to the average weight $`w_{\mathrm{max}}/w`$. This parameter is essential, if we want to transform $`w`$-ted events into $`w=1`$ events, at the latter stage of the MC generation. Minimizing maximum weight is not the same as minimizing variance $`\sigma =\sqrt{w^2w^2}`$. Usually minimizing $`w_{\mathrm{max}}`$ is more difficult. In Foam minimizing variance is also implemented and optionally available. It can be useful in the case when a MC with $`w`$-ted events is acceptable. In fig. 1 we show two examples of the weight distribution evolution in the Foam, when adding more and more cells.

In the Foam algorithm we define two auxiliary distributions $`\rho ^{}(x)`$ and $`\rho _{loss}(x)`$ related to integrand $`\rho (x)`$. Both are constructed together step by step, simultaneously with the construction of the foam of cells, in the exploration process. When our ultimate aim is to minimize $`w_{\mathrm{max}}`$, we define

$$\begin{array}{cc}\hfill \rho ^{}(x)& \underset{yCell_I}{\mathrm{max}}\rho (y),\text{for}xCell_I,\hfill \\ \hfill R_{loss}& =d^nx[\rho ^{}(x)\rho (x)]=d^nx\rho _{loss}(x).\hfill \end{array}$$

here $`\rho _{loss}`$ is the difference between the “ceiling distribution” $`\rho ^{}`$ and the actual distribution $`\rho `$ from which it is derived. The rejection rate in the final MC run is just proportional to $`\rho _{loss}(x)`$ by construction, i.e. the rejection rate $`=R_{loss}/R`$. Also, $`\rho _{loss}(x)`$ has a clear geometrical meaning, see below. When our principal aim is to minimize the ratio of the variance to the average of the weight, $`\sigma /w`$, in the final MC generation, we are led to the following definition:

$$\begin{array}{cc}& \rho ^{}(x)\sqrt{\rho ^2_I},\text{for}xCell_I,\hfill \\ & \rho _{loss}(x)\sqrt{\rho ^2_I}\rho _I,\text{for}xCell_I.\hfill \end{array}$$

The average $`\mathrm{}_I`$ is over the $`I`$-th cell333See ref.  for a detailed derivation of the above prescription.. The ratio $`\sigma /w`$ in the final MC generation is a monotonous ascending function of $`R_{loss}=\rho _{loss}(x)𝑑x^n`$ – hence, minimization of $`R_{loss}`$ is equivalent to minimization of $`\sigma /w`$.

The basic rules governing binary split of a cell is that each split of a cell: $`\omega \omega ^{}+\omega ^{\prime \prime }`$ should decrease the total $`R_{loss}`$ and the decrease should be as big as possible. How to get the best total decrease $`\mathrm{\Delta }R_{loss}`$?
(1) For each next cell split we choose a cell with the biggest $`R_{loss}`$.
(2) Position/direction of a plane dividing a parent cell into two daughter cells is chosen to get the smallest total $`R_{loss}`$.
The actual procedure relies on the small MC exercise done for each single cell during its exploration, in which events are generated with the flat distribution, weighted with $`\rho `$ and projected onto $`n`$ (simplical case) or $`n(n+1)/2`$ (hyperrectangular case) of the cells. Resulting projection histograms are then analysed and the best “division geometry” is found, for which the estimate of $`\mathrm{\Delta }R_{loss}`$ is calculated. How the geometry of the division algorithm looks like? In case of a $`n`$-dimensional simplex defined by $`n+1`$ vertices, a pair of vertices $`x_i`$ and $`x_j`$ is chosen and a new vertex $`Y`$ is put somewhere on the line in between: $`Y=\lambda x_i+(1\lambda )x_j,\mathrm{\hspace{0.33em}\hspace{0.33em}0}<\lambda <1`$. Two daughter simplices are defined with two list of vertices: $`(x_1,\mathrm{},x_{i1},Y,x_{i+1},\mathrm{},x_{j1},x_j,x_{j+1},\mathrm{},x_n,x_{n+1})`$ and $`(x_1,\mathrm{},x_{i1},x_{i+1},\mathrm{},x_{j1},Y,x_{j+1},\mathrm{},x_n,x_{n+1})`$. Now, how do we choose a pair $`(i,j)`$ and the value of $`\lambda `$? A short sample of the MC events (100-1000) is generated $``$ cell. Each MC point is projected (see fig. 2) $`XY`$ onto an edge $`(i,j),ij`$. For each $`(i,j)`$ the $`dN/d\lambda `$ is histogrammed, and its “loss” functional $`R_{loss}`$ is estimated. The optimal edge $`(i,j)`$ with the biggest loss is selected. For the optimal edge the cell division ration $`\lambda `$ is deduced from its histogram ($`\lambda `$ is always a rational number, $`n/N_{bin}`$).

In fig. 3 we show several examples of the evolution of the 2-dimensional foam of cells in case of simplices (triangles) and rectangles for up to 2500 cells.

The important practical question is whether hyperrectangular cells or simplical cells are beter as a building blocks. There is no simple answer. Simplices are limited to low dimensions $`n<6`$ because calculation of increasingly large number of determinants leads to much CPU time consumption. In addition, for simplices memory consumption is $``$ 16$`n`$ Bytes/Cell while for hyperrectangles memory consumption can be limited to below 50Bytes/Cell independently of $`n`$. See ref.  for details. Experimenting with many testing functions has also shown that quite often hyperrectangles provide final MC efficiency as good or sometimes even better than for simplices.

The method of cell division relying on the analysis of the projection histograms requires many function calls in the cell exploration. The overall efficiency of the final MC is best for large number of cells. That points towards necessity of the large number of the function calls during the exploration, which may lead to excessive CPU time consumption. The following efficient trick was found to limit number of the function calls in the cell exploration: During MC exploration of a given new cell we continuously monitor the number of accumulated effective $`W=1`$ events, $`N_{eff}=\frac{(w_i)^2}{{\scriptscriptstyle w_i^2}}`$, and stop exploration when $`N_{eff}/n_{bin}>25`$, where $`n_{bin}`$ is the number of bins in each histogram used to estimate the best division direction/edge and parameter. In this way the increase of $`N_{samp}`$ is not wasted for cells inside which the integrand is already varying very little.

Programs
At present the main development version of the Foam implementation is available in the C++ language. Two other versions written in Fortran 77 are also available. They are functionally equivalent to version in C++444 Future development will be restricted to C++..

Conclusions
Foam algorithm is a versatile adaptive, general-purpose type of an algorithm based on the cellular division of the integration domain. The geometry of the “foam of cells” is rather simple, cells of simplical and/or hyperrectangular shape are constructed in the process of a binary split. It works in principle for arbitrary distribution – no assumption of factorizability as in VEGAS of Ref.. Encoding cells with memory-efficient methods allows up to $`10^6`$ cells to be built in the computer memory of a typical desktop computer.


text:
# Acknowledgements

## Acknowledgements

###

I am very much indebt to Professors I.D. Soares and P.S.Letelier for helpful discussions on the subject of this paper. Special thanks go to Professor Jerry Griffiths for his kind attention and interest in our work. Thanks are also due to Professor Paul Tod for suggesting me this line of work many years ago during my post-doctoral at the Mathematical Institute of Oxford University. Grants from CNPq (Ministry of Science of Brazilian Government) and Universidade do Estado do Rio de Janeiro (UERJ) are acknowledged.


text:
# Status of and performance estimates for QCDOC Presented by PAB at Lattice 2002, Boston.

## Introduction

QCDOC was designed for a cost-effective balance between memory bandwidth, floating point performance and communication performance when running massively parallelised double precision QCD codes. As shall be seen, our simulations suggest QCDOC should scale to exceptionally large machine sizes on a fixed scientifically interesting lattice volume - or equivalently QCDOC will operate efficiently on very small local volumes. This gives the focussed compute power necessary to simulate much lighter quark masses than are accessible today, at a cost of 1 US-dollar per sustained Megaflop.

## Hardware Overview

The QCDOC design is based on an application specific integrated circuit, or ASIC. We use IBM System-On-a-Chip technology to integrate all the node logic on a single silicon chip. The CPU core is a PowerPC 440 whose FPU can execute one operate instruction per clock cycle. Peak performance of two flops per clock is obtained by using fused multiply-add instructions. The on-chip features include 4MByte of embedded DRAM memory and its controller, a memory controller for off-chip DDR memory, an Ethernet interface, and a custom serial communications unit (SCU) linked by a number of on-chip busses.

The prefetching edram controller (PEC) has been custom designed by our colleagues at IBM Research, and provides ample bandwidth for the FPU. It automatically prefetches two read streams, and is multi-ported allowing remote communication to overlap local computation without contention.

The SCU hardware implements a 6 dimensional hyper-torus, with bi-directional low latency links. The aggregate 1.5GByte/s inter-node bandwidth easily supports QCD with 1 Gflop of local floating point all the way down to $`2^4`$ local volumes. The DMA engines in the SCU transfer complex patterns between nodes while the local floating point calculation proceeds uninterrupted.

The high dimensionality of the network both increases the communication bandwidth of the machine, and gives very symmetric local volumes when distributing the problem over very many nodes in four, or even five, dimensions. This yields the most favorable surface to volume ratio.

## Hardware Status

The QCDOC ASIC design is functionally complete, and the preliminary design release is undergoing timing driven physical layout at IBM Raleigh. Verification is on-going with release to manufacture expected by mid fall and first silicon at the end of this year. Large prototype machines are expected in spring 2003.

## Software Environment

Two widespread and standard compilers are being used, namely the IBM xlc compiler and the GNU gcc C and C++ compilers. Debug facilities are provided by a port of the full featured RISCWatch remote debug tool.

The standard runtime support libraries are a popular GNU Public License libc released by Cygnus111Popularised in the Cygwin software package. Runtime libraries are available for performing communication over the physics network via the SciDac QMP interface which is in turn implemented on top of a native SCU interface.

Assembler optimised QCD kernels will be freely available.

## Node kernels

The use of the PowerPC means that standard operating systems could in principle be run on the nodes. However, to scale to very small local volumes it is essential that, in addition to very low latency hardware, we avoid unnecessary software latency.

We will use a lean node kernel to avoid scheduler overhead and use the PowerPC virtual memory (VM) hardware to protect, but not translate, memory pages. This gives both the benefits of robust and graceful error recovery, and the benefits of zero-copy SCU transfers without VM gymnastics.

## Host operating system

The host for QCDOC will be an SMP Unix server with multiple Gigabit links to the Boot/Diagnostic/IO Ethernet network. A multi-threaded qdaemon will boot and manage the machine partitions, and service socket-based connections to these partitions from a number of client programs. Sufficient functionality has already been implemented to boot and download code to PowerPC boards in parallel.

## Performance in simulation

As part of the design verification programme a number of benchmarks and stress tests have been written, to both check functional correctness and that design goals are met. We shall present performance figures for some of these tests, where we have taken a “nominal” 500MHz CPU, and made reasonable assumptions about the interconnect wire length.

Table 1 shows a sample of common vector operations performed in Lattice QCD codes. Most of these figures were obtained using a C++ assembler code generator written by the speaker to automate loop unrolling, prefetching and detailed scheduling.

Table 2 shows performance measurements of optimised implementations of various preconditioned Fermion operators in double precision. In these operators unpaired adds and multiplies set an upper limit (e.g. 780Mflops for the Wilson kernel) somewhat below the “raw” peak. The communication overhead is entirely taken into account, and the performance is excellent despite having to communicate each site over four wires in the $`2^4`$ case.

## Scalability

Performance characteristics at small local volumes are critical to scaling. Figure 1 shows the estimated performance of QCDOC on double precision Wilson dslash code (omits global summation time) distributed in four dimensions. At the smallest volumes the entire operation takes roughly 20 $`\mu `$s, with eight communications in this time. For comparison we overlay a Myrinet Alpha 21264 cluster running single precision assembler code distributed in three dimensions.

Given the fast linear algebra and the scalable matrix multiply the remaining hurdle to be overcome by an iterative solver is global summation. The SCU has “pass-thru” hardware assist for global sums and broadcasts. The global sum has been benchmarked in simulation, and we use the mixture of matrix multiplies, linalg and gsums in the CG algorithm to predict the performance on large machines on Wilson HMC in Table 3.

## Conclusions

QCDOC is progressing well with large machines due in 2003. Simulations indicate these machines will run QCD at high efficiency on the largest machines, while having a low cost per sustained Megaflop and low power consumption.

## Acknowledgements

This research was supported in part by the U.S. Department of Energy, the Institute of Physical and Chemical Research (RIKEN) of Japan, and the U.K. Particle Physics and Astronomy Research Council.


text:
# Measurement of hadronic cross section and preliminary results on the pion form factor using the radiative return at DAΦNE

## 1 Measuring the hadronic cross section with the KLOE detector

### 1.1 Motivation

The measurement of the hadronic cross section at low energy is of great importance for the improvement of the theoretical error of the anomalous magnetic moment of the muon, $`a_\mu =(g_\mu 2)/2`$. The hadronic contribution $`a_\mu ^{hadr}`$ is given by the hadronic vacuum polarization and cannot be calculated at low energy in the framework of perturbative QCD. Following a phenomenological approach, the hadronic contribution can however be evaluated from the measurement of R through a dispersion relation.

$$a_\mu ^{hadr}=(\frac{\alpha m_\mu }{3\pi })^2_{4m_\pi ^2}^{\mathrm{}}𝑑s\frac{R(s)\widehat{K}(s)}{s^2},$$
(1)
where $`R(s)=\frac{\sigma (e^+e^{}hadrons)}{\frac{4\pi \alpha ^2(s)}{3s}}`$ and the kernel $`\widehat{K}(s)`$ is a smooth bounded function growing from $`0.63`$ at threshold to $`1`$ at $`\mathrm{}`$. Due to the $`1/s^2`$ dependence in the integral, hadronic data at low energies are strongly enhanced in the contribution to $`a_\mu ^{hadr}`$. The error of the hadronic contribution is therefore given by the limited knowledge of hadronic cross section data. This error is the dominating contribution to the total error of $`a_\mu ^{theo}`$ ($`\delta a_\mu ^{theo}\delta a_\mu ^{hadr}`$

The discrepancy between the theoretical calculation , and the recent experimental value of $`a_\mu `$  depends if $`\tau `$ data are used or not in the evaluation of $`a_\mu ^{had}`$ ; a recent re-evaluation of the hadronic contribution to $`a_\mu `$  found: $`a_\mu ^{theo}a_\mu ^{exp}`$ = $`(33.9\pm 11.2)\times 10^{10}`$ $`[e^+e^{}based]`$, $`a_\mu ^{theo}a_\mu ^{exp}`$ = $`(16.7\pm 10.7)\times 10^{10}`$ $`[\tau based]`$, corresponding to 3.0 and 1.6 standard deviations respectively, reflecting also a $`1.6\sigma `$ discrepancy between $`e^+e^{}`$ and $`\tau `$ evaluations of $`a_\mu ^{had}`$. A discussion on this subject can also be found in .

This unclear situation makes a precise measurement of the hadronic cross section at low energy mandatory, particularly in the channel $`e^+e^{}\rho \pi ^+\pi ^{}`$, which is the main ingredient for $`a_\mu ^{had}`$. Such a cross section has been already measured by the CMD-2 collaboration with an energy scan in the range below 1.4 GeV ; an accuracy of 0.6% was recently achieved in the region between 0.61 and 0.96 GeV .

### 1.2 Radiative Return

DA$`\mathrm{\Phi }`$NE  is an $`e^+e^{}`$ storage-ring collider working at the $`\varphi `$ resonance (1020 MeV). As an experiment at a collider with a fixed centre of mass energy, KLOE can measure the hadronic cross section $`\sigma `$($`e^+e^{}`$ hadrons) as a function of the hadronic system energy using the radiative return  i. e. studying the process $`e^+e^{}`$ hadrons$`+\gamma `$. The emission of one photon before the beams interact (Initial State Radiation, ISR in the following) lowers the interaction energy and makes possible to produce the hadronic system with an invariant mass varying from the $`\varphi `$ mass down to the production threshold.

The method represents an alternative approach to the conventional energy scan used so far for hadronic cross section measurements. A very solid theoretical understanding of Initial State Radiation (described by the radiation function $`H`$) is mandatory in order to extract the cross section $`\sigma `$($`e^+e^{}`$ hadrons) as a function of $`Q_{had}^2`$ from the measured differential cross section $`d\sigma `$($`e^+e^{}`$ hadrons$`+\gamma `$)$`/dQ_{had}^2`$:

$`Q_{had}^2{\displaystyle \frac{d\sigma (e^+e^{}\mathrm{hadrons}+\gamma )}{dQ_{had}^2}}`$ $`=`$
$`\sigma (e^+e^{}\mathrm{hadrons})H(Q_{had}^2,\theta _\gamma )`$
where $`Q_{had}^2`$ and $`\theta _\gamma `$ are respectively the invariant mass squared of the hadronic system and the acceptance cut on the polar angle of the photon.

Radiative corrections for $`e^+e^{}\pi ^+\pi ^{}\gamma `$, have been calculated up to ISR NLO by different theoretical groups , and recently implemented in a new Monte Carlo generator named PHOKHARA , This generator with a claimed accuracy of 0.5% has been used in the present analysis.

We want to stress that the radiative return has the merit compared with the conventional energy scan that the systematics of the measurement (e.g. normalization, beam energy) are the same for any experimental point and must not be evaluated at each energy step. However a precise determination of the angle of the hard photon as well as the full control of events with the photon emitted in the final state (pure QED FSR or other resonant processes, as for example $`\varphi \pi ^+\pi ^{}\gamma `$ ), is required.

## 2 Analysis of the process $`e^+e^{}\pi ^+\pi ^{}\gamma `$

Due to the importance of the $`2\pi `$ final state for $`a_\mu ^{had}`$, we concentrated at KLOE on the analysis of the radiative process $`e^+e^{}\pi ^+\pi ^{}\gamma `$, whose feasibility has been extensively studied by Monte Carlo .

### 2.1 The KLOE Detector

KLOE  is a typical $`e^+e^{}`$ multiple purpose detector with cylindrical geometry, consisting of a large helium based drift chamber (DC, ), surrounded by an electromagnetic calorimeter (EmC, ) and a superconducting magnet ($`B=0.52`$ T). The detector has been designed for the measurement of $`CP`$ violation in the neutral kaon system, i.e. for precise detection of the decay products of $`K_S`$ and $`K_L`$. These are low momenta charged tracks ($`\pi ^\pm ,\mu ^\pm ,e^\pm `$ with a momentum range from $`150`$ MeV/c to $`270`$ MeV/c) and low energy photons (down to $`20`$ MeV).
The DC dimensions ($`3.3`$ m length, $`2`$ m radius), the drift cell shapes ($`2`$x$`2`$ cm2 cells for the inner $`12`$ layers, $`3`$x$`3`$ cm2 cells for the outer $`46`$ layers) and the choice of the gas mixture ($`90\%`$ Helium, $`10\%`$ Isobutane; $`X_0=900`$ m) had to be optimized for the requirements prevailing at a $`\varphi `$ factory. The KLOE design results in a very good momentum resolution: $`\sigma _p_{}/p_{}0.4\%`$ at high tracking efficiencies ($`>99\%`$).
The EmC is made of a matrix of scintillating fibres embedded in lead, which guarantees a good energy resolution $`\sigma _E/E=5.7\%/\sqrt{E(\mathrm{GeV})}`$ and excellent timing resolution $`\sigma _t=57\mathrm{p}\mathrm{s}/\sqrt{E(\mathrm{GeV})}50`$ ps. The EmC consists of a barrel and two endcaps which are surrounding the cylindrical DC; this gives a hermetic coverage of the solid angle ($`98\%`$). However, the acceptance of the EmC below $`20^{}`$ is reduced due to the presence of quadrupole magnets close to the interaction point and does not allow to measure e.g. the photon of $`\pi ^+\pi ^{}\gamma `$ events with low $`\theta _\gamma `$ angles.

It will be shown in the following that an efficient selection of the $`\pi ^+\pi ^{}\gamma `$ signal is possible, without requiring an explicit photon detection. The relatively simple signature of the signal (2 high momentum tracks from the interaction point, see Fig. 1) and the good momentum resolution of the KLOE tracking detector allow us to perform such a selection.

### 2.2 Signal selection

The selection of $`e^+e^{}\pi ^+\pi ^{}\gamma `$ events is done in the following steps:

* detection of two charged tracks, with polar angle bewteen $`40^o`$ and $`140^o`$, coming from a vertex in the fiducial volume $`R<8`$ cm, $`|z|<15`$ cm. The cuts on the transverse momentum $`p_T>200`$ MeV or on the longitudinal momentum $`|p_z|>90`$ MeV reject tracks spiralizing along the beam line, ensuring good reconstruction conditions. The probability to reconstruct a vertex in the drift chamber is $`95\%`$ and has been studied with Bhabha data, selected using the calorimeter only. The overall tracking reconstruction efficiency has been also evaluated using $`\pi ^+\pi ^{}\pi ^0`$ events selected by detecting $`\pi ^0`$ in the electromagnetic calorimeter;
* identification of pion tracks: a Likelihood Method, using the time of flight of the particle and the shape of the energy deposit in the electromagnetic calorimeter, has been developed to reject the $`e^+e^{}e^+e^{}\gamma `$ background. A control sample of $`\pi ^+\pi ^{}\pi ^0`$ has been used to study the behaviour of pions in the electromagnetic calorimeter, since their interaction are not well reproduced by Monte Carlo. The effect of the selection on $`\pi ^+\pi ^{}\gamma `$ and $`e^+e^{}\gamma `$ events is visible in the upper plot of Fig. 2, as a function of the track mass, $`M_{track}`$: this variable is calculated from the reconstructed momenta, $`\stackrel{}{p}_+`$, $`\stackrel{}{p}_{}`$, applying 4-momentum conservation, under the hypothesis that the final state consists of two particles with the same mass and one photon. For each event class ($`\mu ^+\mu ^{}\gamma `$, $`\pi ^+\pi ^{}\gamma `$) the track mass distribution is peaked at the proper mass; $`\pi ^+\pi ^{}\pi ^0`$ events populate the region on the right of the $`\pi ^+\pi ^{}\gamma `$ peak. As it can be noted from Fig. 2, up, $`e^+e^{}\gamma `$ events are drastically reduced, while the $`\pi ^+\pi ^{}\gamma `$ peak is essentially unaffected by this selection (the estimated signal efficiency is larger than 98%).
* cut on the track mass: $`\mu ^+\mu ^{}\gamma `$ events are rejected by a cut at 120 MeV in the track mass. The discrimination between $`\pi `$ and $`\mu `$ using calorimeter information is not helpful since pions behave frequently like minimum ionizing particles. After this cut we find a contamination of $`\mu ^+\mu ^{}\gamma `$ background smaller than $`1\%`$. $`\pi ^+\pi ^{}\pi ^0`$ events are rejected with a cut in the two-dimensional distribution of the track mass versus the two pion invariant mass squared, $`Q_{\pi \pi }^2`$, shown in Fig. 2, down. Due to the large production cross section for these events ($`BR(\varphi \pi ^+\pi ^{}\pi ^0)15\%`$, $`\sigma 500`$ nb) and the similar kinematics, some residual contamination is expected at small $`Q_{\pi \pi }^2`$ values ($`Q_{\pi \pi }^2<0.4GeV^2`$). The efficiency of the track mass cut, as evaluated from Monte Carlo, is $`90\%`$;
* definition of the angular acceptance: The polar angle of the photon, $`\theta _\gamma `$ is calculated from the charged tracks as $`180^0\theta _{\pi \pi }`$, where $`\theta _{\pi \pi }`$ is the polar angle of the two pion system. Two fiducial volumes for the analysis of $`\pi ^+\pi ^{}\gamma `$ events are defined: $`\theta _\gamma <15^o`$ or $`\theta _\gamma >165^o`$ (small angle analysis) and $`55^o<\theta _\gamma <125^o`$ (large angle analysis). These two regions mainly differ for FSR and background contamination .

In this paper we will show the preliminary results obtained on the analysis at small angle ($`\theta _\gamma <15^o`$ or $`\theta _\gamma >165^o`$, $`40^o<\theta _\pi <140^o`$).

## 3 Final State Radiation

It is now a common understanding that the radiation of one hard photon by final pions should be included in the measurement of the hadronic cross section for the evaluation of $`a_\mu ^{had}`$ .

However , in the case of radiative return, the differential cross section for the process $`e^+e^{}\pi ^+\pi ^{}\gamma `$ is proportional to the pion form factor squared at the measured $`Q_{\pi \pi }^2`$ only for ISR:

$$\left(\frac{\mathrm{d}\sigma }{\mathrm{d}Q^2}\right)_{\mathrm{ISR}}|F_\pi (Q^2)|^2\sigma _{e^+e^{}\pi ^+\pi ^{}}(Q^2),$$
(2)
while for FSR the differential cross section 111If no additional photons are emitted in the initial state. is proportional to the pion form factor evaluated at the $`\varphi `$ mass squared:

$$\left(\frac{\mathrm{d}\sigma }{\mathrm{d}Q^2}\right)_{\mathrm{FSR}}|F_\pi (M_\varphi ^2)|^2\sigma _{e^+e^{}\pi ^+\pi ^{}}(Q^2).$$
(3)
For this reason events due to FSR must be properly rejected by the analysis cuts.

Fig. 3 shows the ratio $`FSR/(ISR+FSR)`$ as a function of the $`Q_{\pi \pi }^2`$ for the small angle region ($`\theta _\gamma <15^o`$ or $`\theta _\gamma >165^o`$, $`40^o<\theta _\pi <140^o`$). As it can be noted the average value of this ratio is below 0.5%.

## 4 Results

### 4.1 Effective cross section

Fig. 4 up shows the number of $`\pi \pi \gamma `$ events selected inside the small angle acceptance cuts at the end of the selection chain, for $`100`$ bins between $`0.02`$ and $`1.02GeV^2`$. This distribution corresponds to $`\mathrm{\hspace{0.17em}73}pb^1`$ of analyzed data out of $`200pb^1`$ collected in 2001 and already reconstructed. More than 1083000 were selected, corresponding to $`15000`$ events/$`pb^1`$.

Even before unfolding the spectrum for the detector resolution effects, the $`\rho \omega `$ interference (as well as the radiative tail) can be clearly seen, showing the excellent momentum resolution of the KLOE DC.

### 4.2 Pion form factor calculation

Neglecting FSR interference, the pion form factor can be extracted as a function of $`Q_{\pi \pi }^2`$, from the observed $`\pi ^+\pi ^{}\gamma `$ spectrum ($`N^{obs}`$) as:

$$|F_\pi (Q_i^2)|^2=\frac{N_i^{obs}}{ϵ(Q_i^2)LH(Q_i^2)}$$
(4)
where $`ϵ(Q_i^2)`$ is the global analysis efficiency, $`L`$ is the integrated luminosity and $`H(Q_i^2)`$ is the NLO cross section for $`e^+e^{}\pi ^+\pi ^{}\gamma `$ (only ISR) under the assumption of pointlike pions.

These quantities were evaluated in the following ways:

1. global efficiency: this is the product of trigger, reconstruction, filtering, likelihood and $`M_{track}`$ cut efficiency. Apart from the latter which was evaluated only by Monte Carlo, the other efficiencies were evaluated from unbiased samples of data with similar kinematics (like $`\pi ^+\pi ^{}\pi ^0`$, Bhabha’s or $`\pi ^+\pi ^{}\gamma `$ itself). The Monte Carlo was used to generate the appropriate $`\pi ^+\pi ^{}\gamma `$ phase space.

Fig. 4, down, shows the dependence of the global efficiency as a function of reconstructed $`Q_{\pi \pi }^2`$ (only statistical errors are shown). Above $`0.4GeV^2`$ the average error is 2%, mainly dominated by the limited Monte Carlo statistics. Below $`0.4GeV^2`$ the error is larger, reflecting again the limited Monte Carlo statistics used in the evaluation of $`M_{track}`$ cut efficiency.
2. Luminosity: The DA$`\mathrm{\Phi }`$NE accelerator does not have luminosity monitors at small angle due to the existence of focusing quadrupole magnets very close to the interaction point. The luminosity is therefore measured using Large Angle Bhabhas ($`55^o<\theta _{+,}<125^o`$, $`\sigma =425`$nb). The number of LAB candidates are counted and normalized to the effective Bhabha cross section obtained from Monte Carlo.

The precision of this measurement depends both on the understanding of experimental efficiencies and acceptances and on the theoretical knowledge of the process.

The systematic errors arising from the LAB selection cuts are well below $`1\%`$. All the selection efficiencies concerning the LAB measurement (Trigger, EmC clusters, DC tracking) are above $`98\%`$ and well reproduced by the detector simulation. The background due to $`\mu ^+\mu ^{}\gamma `$, $`\pi ^+\pi ^{}\gamma `$ and $`\pi ^+\pi ^{}\pi ^0`$ is below $`1\%`$.

KLOE uses two independent Bhabha event generators (the Berends/Kleiss generator  modified for DA$`\mathrm{\Phi }`$NE  and BABAYAGA ).

The very good agreement of the experimental distributions ($`\theta _{+,}`$, $`E_{+,}`$) with the event generators and a cross check with an independent luminosity counter based on $`e^+e^{}\gamma \gamma (\gamma )`$ indicates a precision of better than $`1\%`$.

More systematics checks (e.g. the effect of a varying beam energy and of a displaced beam interaction point) are under way.
3. determination of H: The NLO cross section for pointlike pions $`H`$ was computed for each bin of $`Q_{\pi \pi }^2`$ using a modified version of PHOKHARA, in which $`F_\pi (Q^2)=1`$.

$`210^6`$ events were generated within the acceptance of the small angle analysis: $`\theta _\gamma <15^o`$ or $`\theta _\gamma >165^o`$, $`40^o<\theta _\pi <140^o`$.

Fig. 5 shows a comparison between the differential cross section generated by Monte Carlo with a given parametrization of the pion form factor (white dots), and the same assuming pointlike pions (black dots).

The systematical error on the numerical evaluation of $`H`$ for the given acceptance cuts is estimated to be less than $`1\%`$ by comparing the pion form factor obtained by dividing bin-by-bin the two distributions of Fig. 5 with the expected analytical parametrization.

The pion form factor extracted using eq. 4, is shown in Fig. 6. The spectrum can be divided in three regions according to the measurement error: below $`0.4GeV^2`$ with an error of 5-10%, between $`0.4GeV^2`$ and $`0.5GeV^2`$, with an average error of 3%, and above $`0.5GeV^2`$ with an average error of 2%. The measurement error is dominated by the limited Monte Carlo statistics used in the evaluation of the efficiencies.

### 4.3 Fit of the pion form factor

A preliminary fit was applied to data using the parametrization for $`F_\pi (Q^2)`$ found in  (KS):

$$F_\pi (Q^2)=\frac{BW_\rho \frac{(1+\alpha BW_\omega )}{1+\alpha }+\beta BW_\rho ^{}}{1+\beta }$$
(5)
The mass and the width of the $`\rho `$ as well as $`\alpha `$ and $`\beta `$ were free parameters of the fit, while the other paramaters were kept fixed (to the values of ).

The results obtained from the fit were: $`M_\rho =(772.6\pm 0.5)MeV`$, $`\mathrm{\Gamma }_\rho =(143.7\pm 0.7)MeV`$, $`\alpha =(1.48\pm 0.12)10^3`$, $`\beta =0.147\pm 0.002`$. These values, even if preliminary, are in good agreement with the ones found by CMD-2  using a light different parametrization (Gounaris-Sakurai). This agreement can be also confirmed from Fig. 6, which shows the results of the two parametrization (KLOE blue line, CMD-2 red line) overimposed to our experimental data. The discrepancy at low $`Q_{\pi \pi }^2`$ could be due either to the model inadequacy or to residual background events.

## 5 Conclusion and outlook

The status of $`e^+e^{}\pi ^+\pi ^{}\gamma `$ analysis and the preliminary results obtained on the pion form factor were presented. The analysis was based on an integrated luminosity of $`73pb^1`$ , which corresponds to about 1/3 of the data collected (and already reconstructed) in 2001. The average error for single point above $`0.4GeV^2`$ is $`2\%`$, mainly dominated by the the limited Monte Carlo statistics in the efficiency evaluation. FSR and other background channels (as for example Bhabha’s or $`\pi ^+\pi ^{}\pi ^0`$) are kept below 1% by a likelihood method and kinematical cuts of the analysis. Other sources of systematical errors due to luminosity, radiative corrections and resolution are still under study, and should be included into our errors.

A comparison of the pion form factor with the fitted parametrization obtained by the CMD-2 collaboration, shows a good (preliminary) agreement, even without unfolding the spectrum for the detector resolution effects. Such a comparison, even at an accuracy of 2-3%, has become more and more important in view of the $`1.6\sigma `$ discrepancy between $`e^+e^{}`$ and $`\tau `$ data in the recent evaluation of $`a_\mu ^{had}`$ .

In order to improve the accuracy on $`a_\mu ^{had}`$ a final precision for this measurement below $`1\%`$ is needed. This is very important both in the energy region around the $`\rho `$, where the CMD-2 collaboration has recently reached a 0.6% accuracy as well as below $`0.6GeV`$, which contributes to $`a_\mu `$ with about $`100\times 10^{10}`$ and is known with a worse accuracy . KLOE can study both these regions using $`\pi ^+\pi ^{}\gamma `$ events with the hard photon emitted at small and large angle. However in order to reach such an highly demanding accuracy many factors must be well under control, both experimentally and theoretically. The results here presented, together with the larger statistics up to now collected by KLOE ($`500pb^1`$) and the intense theoretical work on radiative corrections from different groups, are a promising indication for achieving such a challenging task.