text:
# Detection of polarized mm and sub-mm emission from Sgr A⋆

## 1 Introduction

Observations of stellar proper motions in the vicinity of Sgr A (Ghez et al. 1998), the non-thermal radio source at the apparent centre of the Galaxy, show that its mass of $`2.5\times 10^6`$M is highly compact on a scale $`<0.01`$pc, reinforcing its claim to be a massive black-hole candidate. The spectrum of Sgr A extends as a rough power law $`\nu ^\alpha `$, where $`\alpha `$ lies between $`\frac{1}{4}`$ and $`\frac{1}{3}`$, from a low-frequency turn-over at a few GHz up to $``$ 100GHz, (e.g. Mezger, Duschl and Zylka 1996 and references therein), and above this frequency there is evidence of a mm to sub-mm excess over this power law extending almost to the atmospheric cut-off near 1000 GHz (Serabyn et al. 1997, Falcke et al. 1998). The nature of the mm/sub-mm excess has been discussed by Serabyn et al. (1997) and the possibility that it is due to dust is effectively eliminated: at 1.3mm its size is less than 1 arcsecond with an implied brightness temperature in excess of 100 K, inconsistent with the shorter wavelength data, and the implication is that there are two separate synchrotron components.

The radio spectrum of Sgr A has for some years been modelled in various ways in terms of synchrotron radiation but searches for linear polarization, the characteristic signature of this mechanism, have only recently been reported. Bower et al. (1999a,b) used the VLA at 4.8, 8.4, 22 and 43 GHz in its spectropolarimetric mode and the BIMA array at 86 GHz finding upper limits of $`<`$0.2% at the lower frequencies and $`<`$1% at 86 GHz. These small upper limits on linear polarization, given the sensitivity of the observations to large rotation measures due to Faraday rotation (RM $`10^7`$ rad m-2 at 8.4 GHz), are difficult to account for and this is discussed at length by Bower et al. (1999a). Apart from this there has been only a marginal detection at 800 µm of $`4.9\pm 3.2\%`$ at 129 $`\pm 19^{}`$ (Flett and Murray 1991) and a report of circular polarization at 4.8 GHz by Bower et al. (1999c).

In mm to sub-mm polarimetric imaging studies of the central 15 pc of the Galaxy we have detected linear polarization in Sgr A at 750, 850, 1350 and 2000 µm which we present here. Discussion of the field distribution in the circum-nuclear disk (CND) and neighbouring molecular clouds as revealed by these observations will be presented separately.

## 2 Observations and Results

The observations were made using the SCUBA camera (Holland et al. 1999) and its polarimeter (Greaves et al. 1999) on the 15m James Clerk Maxwell Telescope (JCMT) in Hawaii. SCUBA observes simultaneously with arrays of 37 and 91 bolometers, at either 850 and 450 $`\mu `$m respectively, or 750 and 350 $`\mu `$m with a change of filters. Polarimetric observations were made in imaging mode at all four submillimetre wavelengths, and single-point polarimetry towards Sgr A\* was also done at 1.35 and 2mm, using single bolometers and supplementing the data with small grid maps. The arrays have fields of view of 2.3 arcmin, and full-width half-maximum beam sizes range from 7–8′′ at the shortest wavelengths to 34′′ at 2 mm (Table 1).

Observations were made by chopping the secondary mirror at 7.8 Hz to remove the mean sky level, and nodding between left and right beams at slower rates to take out sky gradients. For polarimetry, a rotating half-wave plate and fixed etched grid were used to modulate the signal seen by the detectors. Polarimetry at 750 and 850 $`\mu `$m resulted in a Nyquist-sampled image after 32s of integration, after which the waveplate was stepped by 22.5. The single-pixel data were observed in a similar manner but with 8s of integration per waveplate angle. Total integration times were 5–10 complete waveplate cycles, or 15–85 minutes. The data were reduced using the SURF (Jenness and Lightfoot 1998) and POLPACK (Starlink user note 223) software packages; for each pixel a sinusoidal modulation is fitted to deduce the source percentage and direction ($`p,\theta `$) of the polarization.

From observations of Saturn and Uranus, during the run and from archival data, the instrumental polarization is found to be close to 1% and known to $`\pm 0.1`$$`0.2`$ % at a given wavelength and elevation. At 450 and 1350$`\mu `$m the instrumental polarizations are larger (3.5 and 2.7% respectively due to wind blind structure) but still known to 0.2 to 0.8% respectively. The planets were assumed to be unpolarized and a limited check was made on this for Saturn where sky rotation can be used to separate instrumental from any planetary polarization. While we cannot exclude the possibility that Saturn is slightly polarized this introduces an error of less than 0.3% in the instrumental polarization subtraction. The instrumental polarization changes by only 0.1% over the elevation range used and the uncertainties in the corrected polarization is considered good to $`\pm 0.10.2`$%. Errors in the absolute orientation of the waveplate are one degree or less.

Intensity calibration at 1350 and 2000$`\mu `$m was obtained using Uranus taking the brightness temperatures as 96K and 110K respectively, and at 850$`\mu `$m from Mars taking $`T_B`$ = 208K. Since the emphasis is on the polarimetric results, we did not perform extremely accurate calibration observations, but typical uncertainties are only about $`\pm 10`$ %. For the 750 $`\mu `$m data Saturn was used for calibration taking $`T_B`$=123K and a greater error was introduced since this planet is larger than the beam; we note that the 750 $`\mu `$m flux for Sgr A appears anomalously low. Atmospheric tranmission was measured with skydips.

A fundamental limitation is the maximum chop throw of 180′′ (in this case at a position angle of 145 east of north), which proved to be insufficient at 350 and 450 $`\mu `$m. Surrounding dust emission has a steep spectral index (Pierce-Price et al.. 2000) resulting in increasing off-beam contamination at these shortest wavelengths where Sgr A was not detected above background. At 850 $`\mu `$m the off-beam signals are only 15% of the Sgr A flux and these uncertainties should be similar or smaller at 750, 1350 and 2000 $`\mu `$m at all of which Sgr A was clearly detected. During the observation period the field of view rotates but the chop orientation on the sky is maintained; thus every spatial point in the map retains the same chop positions in $`\alpha `$ and $`\delta `$ although different array bolometers may be involved. The results of polarimetric imaging in the vicinity of Sgr A are presented in Table 1 and the total flux at the position of Sgr A is given in column 3.

There will in general be two contributions to the background, namely from cool dust in the inner cavity and the edges of the CND and free-free emission from the HII region Sgr A West, predominantly from the central ionized filaments.

Dust and free-free emission contribute significantly to the observed flux but in this wavelength region there is no independent evidence at sufficient resolution to determine the contribution of the dust emission component in the central beams. In the sub-mm the flux in annuli from one to two beam radii were determined from the images and at 1350 and 2000$`\mu `$m observations in eight beams circumferentially distributed about SgrA sampled the ambient flux between 1.7–3.7 beam radii from SgrA. Estimates of the free-free contribution in the central beams and these peripheral areas were made from the 2 arcsecond resolution 3.6cm radio continuum maps of Roberts and Goss (1993), after removing the Sgr A point source itself. These fluxes were then scaled to the present wavelengths as $`\nu ^{0.1}`$, appropriate for optically thin free-free emission, and the approximation is made that the central dust contribution is just that derived from the peripheral regions after subtraction of the free-free component (column 7 in the table). The flux from SgrA after subtraction of the dust and free-free contributions (column 8) will also include any local dust excess (or deficit) over background. Since at 2mm the correction for free-free is large due to the large beam size the reliability of the non-thermal flux estimate at this wavelength will be affected.

The observed E vector polarizations and position angles (north through east) are listed in Table 1 in columns 9 and 10. At 450$`\mu `$m there is general polarization over the central 30 arcseconds with no significant change at the position of Sgr A and its average $``$ 3% at a position angle of 100 (E vector) is attributed to dust. At wavelengths between 750 and 2000$`\mu `$m the polarization at the position of Sgr A is well detected at the 10-$`\sigma `$ level. In the sub-mm it stands out from its surroundings and is clearly associated with the flux peak at the position of Sgr A. This can be seen in Fig 1 which shows the central 80 $`\times `$ 100′′ region about Sgr A at 850$`\mu `$m. At 1350 and 2000$`\mu `$m the observations are from single pointed observations where the HII region Sgr A West and Sgr A are the dominant sources in the beam and Sgr A is the only significant contributor to polarization.

## 3 Interpretation

At 2000$`\mu `$m, and to a lesser extent at 1350$`\mu `$m, there are only two contributors to the flux: diffuse free-free emission and Sgr A itself. Since the former is unpolarized the polarization observed confirms the synchrotron nature of Sgr A; it also suggests that the polarization is only associated with the mm/sub-mm excess since it has not been observed at wavelengths of 3.5mm and longer where the power-law spectrum dominates. Note that at 2000$`\mu `$m the largest contributer is from free-free emission and therefore the intrinsic polarization of SgrA must be large.

A large difference of position angle between the mm and sub-mm observations is evident in Table 1 and various ways of understanding this are next investigated.

Although the mm position angles themselves are very close the three shortest wavelength position angles define a straight line when plotted against $`\lambda ^2`$; if a rotation of $`\pi `$ is subtracted from the 2mm angle a rough linear relation is maintained. However the rotation measure then indicated is so large ($`1.44\times 10^6`$ rad/m2 = $`\pi /(\lambda _2^2\lambda _{1.35}^2)`$) that polarization within the band pass of 40GHz at 2mm would be undetectable. Therefore the closeness of the mm position angles indeed indicate that there is little if any Faraday rotation.

Changes in polarization and position angle with time are well documented in the compact cores of blazars and some AGN (e.g. Saikia and Salter 1988) and also have been observed at mm wavelengths (Nartallo et al. 1998). Sgr A is known to be variable on a time scale of months. While there is a separation of some months between the observations presented here, those at 750, 1350 and 2000$`\mu `$m were all made in August 1999 within a week and these do show the large position angle difference between the mm and sub-mm polarizations. Furthermore, the 850$`\mu `$m observation in March also lies on the same trend of position angle with wavelength. Variability is therefore unlikely to have been a significant influence on the position angle, although it seems possible that Sgr A may have increased in brightness during this period.

Since the sub-mm dust emission is polarized it will affect the observed polarization and it is tempting to try to attribute the whole position angle change to dust.

Denoting by $`I_o,I_d`$ and $`I_s`$ the observed central flux, the estimated dust and Sgr A components respectively and by $`q_o,q_d`$ and $`q_s`$ their Stokes $`q`$ parameters, the synchrotron component $`q_s`$ can be found from

$$I_oq_o=I_dq_d+I_sq_s$$

and similarly for the Stokes $`u`$ components.

It is difficult to estimate the dust contribution to polarization from the 750 and 850$`\mu `$m observations because of the proximity of Sgr A and its high intrinsic polarization. Instead we use the 450$`\mu `$m observations in which Sgr A is not detected, and the polarization is $``$ 3% at a position angle of $`100^{}`$. There are some problems here with polarized flux in the reference beams but other indicators of dust alignment in the central parsec are consistent with this estimate: observations at 350$`\mu `$m (Novak et al. 2000) yield 1.9% at 87 near the position of SgrA, and at 100$`\mu `$m the results are similar (Hildebrand et al. 1993) though both these are in a larger beam; the only information at high resolution is from the mid-infrared (Aitken et al. 1996) where integrations over a central 14 arcsecond diameter yield 2.4% at 125, admittedly weighted to warm dust.

The intrinsic polarization resulting from removing the unpolarized free-free contribution and adopting 3% at 100 for dust emission is shown in Table 1 columns 11 and 12. At mm wavelengths these numbers are insensitive to the uncertainties in the adopted dust polarization and even in the submm the range of dust polarizations given in the previous paragraph changes the intrinsic polarization by $`\pm `$3% at most and less than $`\pm `$3. Similarly uncertainties in the relative fractions of the three flux components has little effect on the intrinsic position angle but changing dilution does introduce errors to the intrinsic polarization fraction and at 750 and 2000$`\mu `$m the upper bounds are poorly defined. In Table 1 the errors to the intrinsic polarization are derived by assuming $`\pm `$15% errors to the free-free and $`\pm `$20% errors to the dust contributions except at 2000$`\mu `$m where the adopted error is $`\pm `$40% (these dust contribution errors are derived from the observed variation in the peripheral flux). It is clear that dust polarization cannot explain the abrupt change of position angle between the mm and sub-mm, rather it is increased and it approaches $`\pi `$/2.

We can turn the above argument around to determine what dust polarization is needed in the sub-mm to cause the position angle shift given the intrinsic synchrotron polarization of the mm results, say $`10`$% at 85. With the 850$`\mu `$m fluxes and polarization given in Table 1 the required dust polarization is then 12.0% at 163, larger than any dust emission polarization yet observed. It might be argued that this could be due to unusual conditions close to Sgr A and an additional dust component in the central beam. Any such emission is already included in the Table 1 entry for Sgr A and because the same polarized intensity ($`450`$mJy) has to be supplied from much less dust emission the required polarization would be in excess of 25%.

An assumption made so far is that throughout the wavelength region of these observations the synchrotron radiation remains either optically thin or self-absorbed. Most models of the synchrotron SED in Sgr A consider the mm/sub-mm excess to be self absorbed but it may well be that at the shortest wavelengths it becomes optically thin. In that case the synchrotron position angle shifts through 90 and becomes close to 175 in the sub-mm. The dust corrected sub-mm points are now discrepant by 6 to 14 and this could be attributed to the crudeness of the correction.

Approximating the synchrotron emitting region by a simple slab with constant energy density and magnetic field we find (e.g. Ginzburg and Syrovatskii 1965, 1969, Pacholczyk and Swihart 1967) that for the polarization to change sign near 1mm then $`\nu _{ssa}`$, the synchrotron self absorption frequency, lies in the range 0.6–0.7mm, close to the peak of the mm/sub-mm excess, and the required fields, densities and energies are similar to those in the compact source models of Beckert and Duschl (1997). The polarization transition of this simple model is quite sharp, as is observed, with the polarization quickly approaching saturation on the long wavelength side, changing sign near 1mm and steadily rising at shorter wavelengths.

A variant on this explanation is that self-absorption results in different optical depths with possibly different field distributions being probed as a function of wavelength. In that case, however, one should encounter diminishing polarization towards shorter wavelengths due to the superposition of a range of field orientations, and this is not observed.

## 4 Discussion

Synchrotron radiation from Sgr A has been detected between 750$`\mu `$m and 2000$`\mu `$m with an intrinsic polarization $``$ 10%. This polarization arises from the spectral region of the mm/sub-mm excess. The present observations are compounded by the presence of polarized emission from dust and dilution by free-free emission in these large beams. Although this has led to a number of possible explanations for the large position angle shift between the mm and sub-mm results all but one of these can be eliminated: the closeness of the position angles at 1.35 and 2mm require that there is little or no Faraday rotation, dust polarization can only produce a small shift which is in the wrong sense, and the shift is not due to variability. Only a transition between optically thin synchrotron radiation in the sub-mm and self-absorption in the mm region appears as a plausible explanation of the shift. Future small beam observations will be needed to confirm this result.

A puzzle re-emphasized here is the lack of polarization at 3.5mm (Bower et al.. 1999a,b) and longer wavelengths in view of the large intrinsic polarization of Sgr A at 1.35 and 2mm. Sufficient dilution of the mm/sub-mm excess to bring it below 1% requires not only a steeply falling spectrum with decreasing frequency for the excess, where $`\nu ^{5/2}`$ is expected for self absorption, but also that the cut-off frequency of the power-law section is below 2mm to ensure a fast rising flux of the power law region. It is possible, as well, that a flare in Sgr A at 100GHz in March 1998 (Tsuboi, Miyazaki and Tsutsumi 1998) has affected the result of Bower et al. 1999b.

The synchrotron self-absorption transition in the sub-mm proposed here is consistent with the model of the excess by Beckert and Duschl (1997), and requires the excess to arise from a compact source $``$ 2 Schwarzschild radii in size. Also because of the transition the $`\pi `$/2 ambiguity of the relationship between position angle and magnetic field is removed and hence the field orientation is close to east-west, since at mm wavelengths the emission is considered self-absorbed.

## 5 Conclusions

We report the observation of millimetre and sub-millimetre polarization from Sgr A, confirming the role of synchrotron radiation. The polarization is a property of the mm/sub-mm excess, demonstrating that the excess is real and not an artefact of variability and that it and the power-law spectrum arise in distinct structures. There is a large position angle shift between the mm and sub-mm observations. This can be explained by a transition between optically thin and self absorbed synchrotron radiation near 1 mm. Such a high self-absorption frequency implies a very compact source $``$ 2RS.

We are grateful to the UK Panel for Allocation of Telescope Time for the award of observing time for this project. JSR acknowledges a Royal Society Fellowship and DP-P a PPARC Research Studentship. The JCMT is operated by the Joint Astronomy Centre on behalf of PPARC of the UK, the Netherlands OSR, and NRC Canada. The authors would like to thank the NCSA Astronomy Digital Image Library (ADIL) for providing images for this paper. We thank an unknown referee for useful comments.


format: text

text:
# dcs: An Implementation of DATALOG with Constraints

## General Info

$`\mathrm{DC}`$ is an answer set programming (ASP) system (?) similar to propositional logic but extended to include Horn clauses. The semantics of $`\mathrm{DC}`$ is a natural extension of the semantics of propositional logic. The $`\mathrm{DC}`$ system is implemented in two modules, ground and dcs, which are written in the ’C’ programming language and compiled with gcc. There are approximately 2500 lines of code for ground and approximately 1500 for dcs. $`\mathrm{DC}`$ has been implemented on both SUN SPARC and a PC running linux.

A $`\mathrm{DC}`$ theory (or program) consists of constraints and Horn rules (DATALOG program). This fact motivates our choice of terminology — DATALOG with constraints. We start a discussion of $`\mathrm{DC}`$ with the propositional case. Our language is determined by a set of atoms $`At`$. We will assume that $`At`$ is of the form $`At=At_CAt_H`$, where $`At_C`$ and $`At_H`$ are disjoint.

A $`\mathrm{DC}`$ theory (or program) is a triple $`T_{dc}=(T_C,T_H,T_{PC})`$, where
1. $`T_C`$ is a set of propositional clauses $`\neg a_1\mathrm{}\neg a_mb_1\mathrm{}b_n`$ such that all $`a_i`$ and $`b_j`$ are from $`At_C`$,
2. $`T_H`$ is a set of Horn rules $`a_1\mathrm{}a_mb`$ such that $`bAt_H`$ and all $`a_i`$ are from $`At`$,
3. $`T_{PC}`$ is a set of clauses over $`At`$.

By $`At(T_{dc})`$, $`At_C(T_{dc})`$ and $`At_{PC}(T_{dc})`$ we denote the set of atoms from $`At`$, $`At_C`$ and $`At_{PC}`$, respectively, that actually appear in $`T_{dc}`$.

With a $`\mathrm{DC}`$ theory $`T_{dc}=(T_C,T_H,T_{PC})`$ we associate a family of subsets of $`At_C(T_{dc})`$. We say that a set $`MAt_C(T_{dc})`$ satisfies $`T_{dc}`$ (is an answer set of $`T_{dc}`$) if
1. $`M`$ satisfies all the clauses in $`T_C`$, and
2. the closure of $`M`$ under the Horn rules in $`T_H`$, $`M^c=LM(T_HM)`$ satisfies all clauses in $`T_{PC}`$ ($`LM(P)`$ denotes the least model of a Horn program $`P`$).

Intuitively, the collection of clauses in $`T_C`$ can be thought of as a representation of the constraints of the problem, Horn rules in $`T_H`$ can be viewed as a mechanism to compute closures of sets of atoms satisfying the constraints in $`T_C`$, and the clauses in $`T_{PC}`$ can be regarded as constraints on closed sets (we refer to them as post-constraints). A set of atoms $`MAt_C(T_{dc})`$ is a model if it (propositionally) satisfies the constraints in $`T_C`$ and if its closure (propositionally) satisfies the post-constraints in $`T_{PC}`$. Thus, the semantics of $`\mathrm{DC}`$ retains much of the simplicity of the semantics of propositional logic. $`\mathrm{DC}`$ was introduced by the authors in (?).

## Applicability of the System

$`\mathrm{DC}`$ can be used as a computational tool to solve search problems. We define a search problem $`\mathrm{\Pi }`$ to be determined by a set of finite instances, $`D_\mathrm{\Pi }`$, such that for each instance $`ID_\mathrm{\Pi }`$, there is a finite set $`S_\mathrm{\Pi }(I)`$ of all solutions to $`\mathrm{\Pi }`$ for the instance $`I`$. For example, the problem of finding a hamilton cycle in a graph is a search problem: graphs are instances and for each graph, its hamilton cycles (sets of their edges) are solutions.

A $`\mathrm{DC}`$ theory $`T_{dc}=(T_C,T_H,T_{PC})`$ solves a search problem $`\mathrm{\Pi }`$ if solutions to $`\mathrm{\Pi }`$ can be computed (in polynomial time) from answer sets to $`T_{dc}`$. Propositional logic and stable logic programming (??) are used as problem solving formalisms following the same general paradigm.

We will illustrate the use of DC to solve search problems by discussing the problem of finding a hamilton cycle in a directed graph. Consider a directed graph $`G`$ with the vertex set $`V`$ and the edge set $`E`$. Consider a set of atoms $`\{hc(a,b):(a,b)E\}`$. An intuitive interpretation of an atom $`hc(a,b)`$ is that the edge $`(a,b)`$ is in a hamilton cycle. Include in $`T_C`$ all clauses of the form $`\neg hc(b,a)\neg hc(c,a)`$, where $`a,b,cV`$, $`bc`$ and $`(b,a),(c,a)E`$. In addition, include in $`T_C`$ all clauses of the form $`\neg hc(a,b)\neg hc(a,c)`$, where $`a,b,cV`$, $`bc`$ and $`(a,b),(a,c)E`$. Clearly, the set of propositional variables of the form $`\{hc(a,b):(a,b)F\}`$, where $`FE`$, satisfies all clauses in $`T_C`$ if and only if no two different edges in $`F`$ end in the same vertex and no two different edges in $`F`$ start in the same vertex. In other words, $`F`$ spans a collection of paths and cycles in $`G`$.

To guarantee that the edges in $`F`$ define a hamilton cycle, we must enforce that all vertices of $`G`$ are reached by means of the edges in $`F`$ if we start in some (arbitrarily chosen) vertex of $`G`$. This can be accomplished by means of a simple Horn program. Let us choose a vertex, say $`s`$, in $`G`$. Include in $`T_H`$ the Horn rules $`hc(s,t)vstd(t)`$, for every edge $`(s,t)`$ in $`G`$. In addition, include in $`T_H`$ Horn rules $`vstd(t),hc(t,u)vstd(u)`$, for every edge $`(t,u)`$ of $`G`$ not starting in $`s`$. Clearly, the least model of $`FT_H`$, where $`F`$ is a subset of $`E`$, contains precisely these variables of the form $`vstd(t)`$ for which $`t`$ is reachable from $`s`$ by a nonempty path spanned by the edges in $`F`$. Thus, $`F`$ is the set of edges of a hamilton cycle of $`G`$ if and only if the least model of $`FT_H`$, contains variable $`vstd(t)`$ for every vertex $`t`$ of $`G`$. Let us define $`T_{PC}=\{vstd(t):tV\}`$ and $`T_{ham}(G)=(T_C,T_H,T_{PC})`$. It follows that hamilton cycles of $`G`$ can be reconstructed (in linear time) from answer sets to the $`\mathrm{DC}`$ theory $`T_{ham}(G)`$. In other words, to find a hamilton cycle in $`G`$, it is enough to find an answer set for $`T_{ham}(G)`$.

This example illustrates the simplicity of the semantics of $`\mathrm{DC}`$ — it is only a slight adaptation of the semantics of propositional logic to the case when in addition to propositional clauses we also have Horn rules in theories. It also illustrates the power of $`\mathrm{DC}`$ to generate concise encodings. All known propositional encodings of the hamilton-cycle problem require that additional variables are introduced to “count” how far from the starting vertex an edge is located. Consequently, propositional encodings are much larger and lead to inefficient solutions to the problem.

## Description of the System

In this section we will discuss general features of $`\mathrm{DC}`$. First, we will discuss the language for encoding problems and give an example by showing the encoding of the hamiltonicity problem. Second, we will describe how we execute the $`\mathrm{DC}`$ system. Third, we will give some details concerning the implementation of dcs, our solver. Last we discuss the expressitivity of $`\mathrm{DC}`$.

In the previous section we described the logic of $`\mathrm{DC}`$ in the propositional case. Our definitions can be extended to the predicate case (without function symbols). Each constraint is treated as an abbreviation of a set of its ground substitutions. This set is determined by the set of constants appearing in the theory and by additional conditions associated with the constraint ( we illustrate this idea later in this section). When constructing predicate $`DC`$-based solutions to a problem $`\mathrm{\Pi }`$ we separate the representation of an instance to $`\mathrm{\Pi }`$ from the constraints that define $`\mathrm{\Pi }`$. The representation of an instance of $`\mathrm{\Pi }`$ is a collection of facts or the extensional database (EDB). The constraints and rules that define $`\mathrm{\Pi }`$ will be referred to as the intensional database (IDB) and the language for writing the problem descriptions that constitute the IDB will be referred to as $`L_{dc}`$. The separation of IDB and EDB means only one predicate description of $`\mathrm{\Pi }`$ is needed.

The modules of $`\mathrm{DC}`$, ground and dcs provide a complete system for describing and finding solutions to problems. An IDB in $`L_{dc}`$ along with a specific EDB are the input to ground. A grounded propositional theory $`T_{dc}`$ is output by ground and used as input to dcs.

A problem description in $`L_{dc}`$ defines predicates, declares variables, and provides a description of the problem using rules. The predicates are defined using types from the EDB. Similarly, the variables are declared using types from the EDB. The rules consist of constraints, Horn rules and post-constraints. The constraints and post-constraints use several constructs to allow a more natural modeling. These constructs could be directly translated to clauses. (We use them as shorthands to ensure the conciseness of encodings.)

We present here a brief discussion of the constraints, Horn rules and post-constraint in $`L_{dc}`$ and their meanings. Let $`PRED`$ be the set of predicates occurring in the IDB. For each variable $`X`$ declared in the IDB the range $`R(X)`$ of $`X`$ is determined by the EDB.

where $`n,m`$ are nonnegative integers such that $`nm,qPRED`$ and $`p_1,\mathrm{},p_i`$ are EDB predicates or logical conditions (logical conditions can be comparisons of arithmetic expressions or string comparisons). The interpretation of this constraint is as follows: for every $`\stackrel{}{x}R(\stackrel{}{X})`$ at least $`n`$ atoms and at most $`m`$ atoms in the set $`\{q(\stackrel{}{x},\stackrel{}{y}):\stackrel{}{y}R(\stackrel{}{Y})\}`$ are true.

where $`n,m`$ are nonnegative integers such that $`nm,qPRED`$. The interpretation of this constraint is as follows: for every $`\stackrel{}{x}R(\stackrel{}{X})`$ at least $`n`$ atoms and at most $`m`$ atoms in the set $`\{q(\stackrel{}{x},\stackrel{}{y}):\stackrel{}{y}R(\stackrel{}{Y})\}`$ are true.

where $`n,m`$ are nonnegative integers such that $`nm,q_1,\mathrm{},q_jPRED`$. The interpretation of this constraint is as follows: for every $`\stackrel{}{x}R(\stackrel{}{X})`$ at least $`n`$ atoms and at most $`m`$ atoms in the set $`\{q_1(\stackrel{}{x}),\mathrm{},q_j(\stackrel{}{x})\}`$ are true.

Certain choices of $`n,m`$ in any of the Select constraints allow construction of even more specific constraints in $`T_{dc}`$:
$`n=m`$ exactly $`n`$ atoms must be true.
$`n=0`$ at most $`m`$ atoms must be true.
$`m=999`$ at least $`n`$ atoms must be true.

where $`q_1,\mathrm{},q_iPRED`$. For every $`\stackrel{}{x}R(\stackrel{}{X})`$ at least one atom in the set $`\{q_1(\stackrel{}{x}),\mathrm{},q_i(\stackrel{}{x})\}`$ must be false.

where $`q_1,\mathrm{},q_iPRED`$. For every $`\stackrel{}{x}R(\stackrel{}{X})`$ at least one atom in the set $`\{q_1(\stackrel{}{x}),\mathrm{},q_i(\stackrel{}{x})\}`$ must be true.

where $`p_1,\mathrm{},p_i,q_1,\mathrm{},q_jPRED`$. For every $`\stackrel{}{x}R(\stackrel{}{X})`$ if all atoms in the set $`\{p_1(\stackrel{}{x}),\mathrm{},p_i(\stackrel{}{x})\}`$ are true then at least one atom in the set $`\{q_1(\stackrel{}{x}),\mathrm{},q_j(\stackrel{}{x})\}`$ must be true. (This constraint represents standard propositional logic clauses.)

where $`p_1,\mathrm{},p_i,q_1,\mathrm{},q_jPRED`$. For every $`\stackrel{}{x}R(\stackrel{}{X})`$ if all atoms in the set $`\{p_1(\stackrel{}{x}),\mathrm{},p_i(\stackrel{}{x})\}`$ are true then all atoms in the set $`\{q_1(\stackrel{}{x}),\mathrm{},q_j(\stackrel{}{x})\}`$ must be true.

where $`p_1,\mathrm{},p_i,q_1,\mathrm{},q_jPRED`$.

### Methodology

Here we show the encoding of a problem in $`\mathrm{DC}`$. We will use the example of hamiltonicity of a graph which we discussed previously. Figure 1 shows the EDB format used by ground. This format is compatible to that used by smodels and others. However, we require that the EDB be in separate files from the IDB. The format for the EDB allows data to be entered as sets, ranges, or individual elements and constant values can be entered on the command line.

The IDB provides a definition of the problem in $`L_{dc}`$. The IDB file has three parts. First a definition of the predicates, next the declaration of variables and last a set of constraints and Horn clauses. The types used in the IDB must be in the data file(s). For example, the only data types in the graph file (Fig. 1) are vtx and edge. Thus the only data types which can be used in the IDB are vtx and edge.

In the idbpred section of Fig. 2, we define two predicates, the vstd predicate and the hc predicate. The vstd predicate has one parameter of type vtx and hc has two parameters both of type vtx.

The idbvar section of Fig. 2 declares two variables $`X`$, $`Y`$ both of type vtx.

The section containing the constraints, Horn clauses, and post-constraints is proceeded by the keyword idbrules (see Fig.2). The order in which the rules are entered is not important. The first constraint, $`Select(1,1,Y;edge(X,Y))hc(X,Y).`$, ensures that each vertex has exactly one outgoing edge. The second constraint, $`Select(1,1,X;edge(X,Y))hc(X,Y).`$, requires that each vertex has exactly one incoming edge.

The first Horn rule ranges over all $`X,Y`$ such that $`edge(X,Y)EDB`$ and $`Xi`$ where $`i`$ is a constant used to initialize $`vstd(i)`$. This rule requires both $`vstd(X)`$ and $`hc(X,Y)`$ to be true before we can assign the value true to $`vstd(Y)`$. The second Horn rule only requires $`hc(X,Y)`$ to be true before $`vstd(Y)`$ is assigned value true; however, in this rule $`X`$ is restricted to the value of $`i`$.

The last line is a post-constraint that requires $`vstd(X)`$ to be true for all $`XR(X)`$ ensuring that the cycle is closed.

### Running the system

Here we will describe the steps for execution of the $`\mathrm{DC}`$ system. The first module of $`\mathrm{DC}`$, ground, has as input a data file(s), a rule file and command line arguments. Output is the theory file which will be used as input to dcs. The theory is written to a file whose name is the catenation of the constants and file names given as command line arguments. The extension .tdc is then appended to the output file name. The command line arguments for ground:

ground -r rf -d df \[-c label=v\] \[-V\]

Required arguments

* rf is the file describing the problem. There must be exactly one rule file.
* df must be one or more files containing data that will be used to instantiate the theory. It is often convenient to use more than one file for data.

Optional arguments

* This option allows use of constants in both data and rule files. When label is found while reading either file it is replaced by v. v can be any string that is valid for the data type. If label is to be used in a range then v must be an integer. For example, if the data file contains the entry queens\[1..q\]. then we can define the constant q with the option -c q=8. If more than one constant is needed then -c b=3 n=14 defines both constants b,n using the -c option.
* The verbose options sends output to stdout during the execution of ground. This output may be useful for debugging of the data or rule files.

For the example of hamiltonicity we could have a data file (see Fig. 1) named 1.gph and an IDB file (see Fig. 2) named hcp. The constant i is needed in the IDB file to initialize the first vertex in the graph. The command line argument would be:

ground -r hcp -d 1.gph -c i=1

The theory file produced would be named 1\_1.gph\_hcp.tdc.

The second module of the $`\mathrm{DC}`$ system, dcs, has as input the theory file produced by ground. A file named dcs.stat is created or appended with statistics concerning the results of executing dcs on the theory. The command line arguments are:

dcs -f filename \[-A\] \[-P\] \[-C\] \[-V\]

Required arguments

* filename is the name of the file containing a theory produced by ground.

Optional arguments

* Prints the positive atoms for solved theories in readable form.
* Prints the input theory and then exits.
* Counts the number of solutions. This information is recorded in the statistics file.
* Prints information during execution (branching, backtracking, etc). Useful for debugging.

### Discussion of dcs

The $`\mathrm{DC}`$ solver uses a Davis-Putnam type approach, with backtracking, propagation and LookAHead to deal with constraints represented as clauses, select constraints and Horn rules, and to search for answer sets. The LookAHead in $`\mathrm{DC}`$ is similar to local processing performed in csat (?). However, we use different methods to determine how many literals to consider in the LookAHead phase. Other techniques, especially propagation and search heuristics, were designed specifically for the case of $`\mathrm{DC}`$ as they must take into account the presence of Horn rules in programs.

Propagation consists of methods to reduce the theory. Literals which appear in the heads of Horn rules, $`l_hAt_H`$ require different interpretations. A literal is a $`l_h`$ if and only if it appears in the head of a Horn rule. We can not guess an assignment for $`l_h`$ rather it must be computed. We can only assign value true to $`l_h`$ if it appears in the head of a Horn rule for which all literals in the body of the Horn rule have been assigned the value true. If one or more literals in the body of a Horn rule have been assigned the value false then that rule is removed. If $`l_h`$ has not been assigned a value and all Horn rules in which $`l_h`$ appeared in the head have been removed then $`l_h`$ is assigned the value false. If we have a post-constraint that required a value be assigned to $`l_h`$ and the value cannot be computed then we must backtrack.

Non Horn rules are constraints which must be satisfied. These rules are identified by tags which indicate which method is needed to evaluate the constraint during propagation.

The LookAHead procedure tests a number of literals not yet assigned a value. The LookAHead procedure is similar to the local processing procedure used for csat (?). A literal is chosen, assigned the value true, then false using propagation to reduce and evaluate the resulting theories. If both evaluations of assignments result in conflicts then we return false and backtracking will result. If only one evaluation results in conflict then we can assign the literal the opposite value and continue the LookAHead procedure. If neither evaluation results in conflict we cannot make any assignments, but we save information (number of forced literals and satisfied constraints) computed during the evaluations of the literal.

The number of literals to be tested has been determined empirically. It is obvious that if all unassigned literals were tested during each LookAHead it would greatly increase the time. However, if only a small number of literals are to be tested during each LookAHead then they must be chosen to provide the best chance of reducing the theory. To choose the literals with the best chance of reducing the theory, we order the unassigned literals based on a sum computed by totaling the weights of the unsatisfied constraints in which they appear. The constraint weight is based on both the length and type of constraint. The shorter the constraint the larger the weight and when literals are removed the weight is recomputed. The literals are tested in descending order of the sum of constraints weights. Using this method we need to test only a very small number of literals during each LookAHead to obtain good results.

At the completion of the LookAHead procedure, we use the information computed during the evaluation of literals to choose the next branching literal and its initial truth assignment.

### Expressitivity

The expressive power of $`\mathrm{DC}`$ is the same as that of logic programming with the stable-model semantics. The following theorem is presented in (?)

###### Theorem 1

The expressive power of $`\mathrm{DC}`$ is the same as that of stable logic programming. In particular, a decision problem $`\mathrm{\Pi }`$ can be solved uniformly in $`\mathrm{DC}`$ if and only if $`\mathrm{\Pi }`$ is in the class NP.

## Evaluating the System

The $`\mathrm{DC}`$ system provides a language, $`L_{dc}`$, which facilitates writing problem descriptions. Once an IDB is written in $`L_{dc}`$ it can be used for any instance of the problem for which data, in the EDB format, is available or can be generated. It is possible to add constraints to IDB for a given problem when new requirements or constraints occur. The constructs used in $`L_{dc}`$ allow for a natural description of constraints. Users need only know enough about a specific problem to be able to describe the problem in $`L_{dc}`$ (there is a user’s manual with examples). To help with programming in $`\mathrm{DC}`$ ground provides error messages and compiling information that are useful for debugging the IDB.

### Benchmarks

The $`\mathrm{DC}`$ system has been executed using problems from NP, combinatorics, and planning. In particular, it has been used to compute hamilton cycles and colorings in graphs, to solve the $`N`$-queens problem, to prove that the pigeonhole problem has no solution if the number of pigeons exceeds the number of holes, and to compute Schur numbers.

The instances for computing hamilton cycles were obtained by randomly generating one thousand directed graphs with $`v=30,40\mathrm{},80`$ and density such that $`50\%`$ of the graphs contained hamilton cycles.

The graphs for instances of coloring were one hundred randomly generated graphs for $`v=50,100,\mathrm{},300`$ with density such that $`50\%`$ had solutions when encode as 3-coloring.

### Comparison

We have used the benchmarks to compare $`\mathrm{DC}`$ with systems based on stable model semantics and satisfiability. The performance of $`\mathrm{DC}`$ solver dcs was compared with smodels, a system for computing stable models of logic programs (?), and csat, a systems for testing propositional satisfiability (?). In the case of smodels we used version 2.24 in conjunction with the grounder lparse, version 0.99.41 (?). The extended rules (?) allowed in the newer versions of smodels and lparse were used where applicable. The programs were all executed on a Sun SPARC Station 20. For each test we report the cpu user times for processing the corresponding propositional program or theory. We tested all three system using the benchmarks discussed in the previous section.

Note that csat performs comparable to $`\mathrm{DC}`$ for pigeonhole (see Fig. 6), N-queens (see Fig. 5) and coloring (see Fig. 4). These are problems where $`\mathrm{DC}`$ encodings do not use Horn rules (in the examples here only encodings for computing hamilton cycles use Horn rules). The closure properties of Horn rules allow for much smaller theories as shown in Table 2. The satisfiability theories for computing hamilton cycles are so large that they were not practical to execute for over $`40`$ vertices (see Fig. 3). smodels performs much better than satisfiability solvers for computing the hamilton cycles although not as well as $`\mathrm{DC}`$. The results for computing Schur numbers (see Table 1) also show much better results for $`\mathrm{DC}`$.

Experimental results show that dcs often outperforms systems based on satisfiability as well as systems based on non-monotonic logics, and that it constitutes a viable approach to solving problems in AI, constraint satisfaction and combinatorial optimization. We believe that our focus on short encodings (see Fig. 2) is the key to the success of dcs.


text:
# A new approach to Ginsparg-Wilson fermions

## Abstract

We expand the most general lattice Dirac operator $`D`$ in a basis of simple operators. The Ginsparg-Wilson equation turns into a system of coupled quadratic equations for the expansion coefficients. Our expansion of $`D`$ allows for a natural cutoff and the remaining quadratic equations can be solved numerically. The procedure allows to find Dirac operators which obey the Ginpsparg-Wilson equation with arbitrary precision.

In the last two years we have witnessed exciting developments for chiral fermions on the lattice (see  for some reviews). It was found that the perfect action approach  leads to the fixed point Dirac operator  which obeys the Ginsparg-Wilson equation . Also the overlap formalism for chiral fermions on the lattice  has led to an explicit expression for a Dirac operator which solves the Ginsparg-Wilson equation, the so-called overlap operator . It was realized, that based on the Ginsparg-Wilson relation it is possible to redefine the chiral symmetry for lattice fermions and this new construction has spurred a wealth of interesting work .

Unfortunately the numerical implementation of the two known solutions of the Ginsparg-Wilson equation seems to be rather troublesome. Computing the inverse square root which appears in the overlap operator is a numerically challenging and expensive task. For the perfect action the situation is even worse, since already the blocking procedure necessary for its computation is extremely costly and so far has been implemented only in two dimensions .

In this letter we present a new approach to the Ginsparg-Wilson equation: The first step is to expand the most general Dirac operator $`D`$ in a basis of simple operators on the lattice. In the second step we insert the expanded $`D`$ into the Ginsparg-Wilson equation which then turns into a system of coupled quadratic equations for the expansion coefficients of $`D`$. Solving the Ginsparg-Wilson equation is equivalent to finding solutions for this system of quadratic equations. Any solution of the Ginsparg-Wilson equation, such as the overlap operator or the perfect action, corresponds to a solution of our system of quadratic equations. Many interesting aspects can be addressed in this new framework, such as e.g. the question if there is a continuous manifold of solutions in the space of coefficients of $`D`$.

Here we concentrate on exploring the possibilities for solving our system of quadratic equations numerically and using this approach for constructing new solutions of the Ginsparg-Wilson equation. We show that our expansion of $`D`$ provides for a natural cutoff which turns the quadratic equations into a simple finite system. Solving this system numerically gives the values for the expansion coefficients. This procedure allows to construct approximate solutions which obey the Ginsparg-Wilson equation with arbitrary precision.

Before we outline the details of our construction, let us begin with some remarks: Usually the derivative term on the lattice is discretized by the following nearest neighbor term ($`U_\mu (x)`$ SU(N) and we set the lattice spacing to 1):

$$\frac{1}{2}\underset{\mu =1}{\overset{4}{}}\gamma _\mu \left[U_\mu (x)\delta _{x+\widehat{\mu },y}U_\mu (x\widehat{\mu })^1\delta _{x\widehat{\mu },y}\right].$$
(1)
However, it is perfectly compatible with all the symmetries to instead discretize the derivative term using

$`{\displaystyle \frac{1}{4}}{\displaystyle \underset{\mu =1}{\overset{4}{}}}\gamma _\mu [U_\mu (x)U_\mu (x+\widehat{\mu })\delta _{x+2\widehat{\mu },y}`$ (2)
$`U_\mu (x\widehat{\mu })^1U_\mu (x2\widehat{\mu })^1\delta _{x2\widehat{\mu },y}],`$ (3)
and there are many more terms that are eligible. Thus an ansatz for the most general Dirac operator $`D`$ must allow for a superposition of all of the possible discretizations for the derivative term. In order to remove the doublers we also have to allow for terms that come with 1I in spinor space and again we allow for all possible terms. We generalize our $`D`$ further, by including all terms also for the remaining elements $`\mathrm{\Gamma }_\alpha `$ of the Clifford algebra, i.e. we include also tensor, pseudovector and pseudoscalar terms. From our examples (1), (3) it can be seen that the terms are characterized by the orientations of the link variables $`U_\mu (x)`$. The ensemble of links supporting the gauge transporters can be viewed as paths and we can denote our ansatz for the most general lattice Dirac operator in the following form:

$$D_{x,y}=\underset{\alpha =1}{\overset{16}{}}\mathrm{\Gamma }_\alpha \underset{p𝒫_{x,y}^\alpha }{}c_p^\alpha \underset{lp}{}U_l.$$
(4)
To each generator $`\mathrm{\Gamma }_\alpha `$ of the Clifford algebra and to each pair of points $`x,y`$ on the lattice we assign a set $`𝒫_{x,y}^\alpha `$ of paths $`p`$ connecting the two points. Each path is weighted with some complex weight $`c_p^\alpha `$ and the construction is made gauge invariant by including the ordered product of the gauge transporters $`U_l`$ $``$ SU(N) for all links $`l`$ of $`p`$. The action is then given by $`S=_{x,y}\overline{\psi }_xD_{x,y}\psi _y`$, where $`x`$ and $`y`$ run over all of the lattice.

The next step is to impose on $`S`$ the symmetries which we want to maintain: Translation and rotation invariance and invariance under C and P. In addition we require our $`D`$ to be $`\gamma _5`$-hermitian, i.e. we require $`\gamma _5D\gamma _5=D^{}`$. This property can be seen to correspond to what leads to the CTP theorem in Minkowski space, i.e. the vector generators $`\gamma _\mu `$ come with a derivative term etc. (compare Equation (13) below).

Translation invariance requires the sets $`𝒫_{x,y}^\alpha `$ of paths and their coefficients to be invariant under simultaneous shifts of $`x`$ and $`y`$, while rotation invariance implies that a term and its rotated image have the same weight. Parity implies that for each path $`p`$ (with coefficient $`c_p^\alpha `$) we must include the parity-reflected copy with coefficient $`s_{parity}^\alpha c_p^\alpha `$ where the signs $`s_{parity}^\alpha `$ are defined by $`\gamma _4\mathrm{\Gamma }_\alpha \gamma _4=s_{parity}^\alpha \mathrm{\Gamma }_\alpha `$.

More interesting are the symmetries C and $`\gamma _5`$-hermiticity. It is easy to see that both of them imply a relation between the coefficient for a path $`p`$ and the coefficient of the inverse path $`p^1`$. Implementing both these symmetries firstly restricts all coefficients $`c_p^\alpha `$ to be either real or purely imaginary. Secondly we find that the coefficient for a path $`p`$ and the coefficient for its inverse $`p^1`$ are equal up to a sign $`s_{charge}^\alpha `$ defined by $`C\mathrm{\Gamma }_\alpha C=s_{charge}^\alpha \mathrm{\Gamma }_\alpha ^T`$, where $`T`$ denotes transposition and $`C`$ is the charge conjugation matrix. To be specific, we use the chiral representation for the $`\gamma _\mu `$ and $`C=i\gamma _2\gamma _4`$. When implementing all these symmetries we find that paths in our ansatz become grouped together where, up to sign factors, all paths in a group come with the same coefficient.

Before we present the explicit form of the most general Dirac operator with the above symmetries, let us first introduce an abbreviated notation: The terms in the Dirac operator are determined by the corresponding element of the Clifford algebra and the paths with their relative signs. For the two examples (1), (3) the paths are given by one (two) steps in positive $`\mu `$-direction and one (two) steps in negative $`\mu `$-direction with a relative minus sign. We now denote an arbitrary path of length $`n`$ on the lattice by $`<l_1,l_2,\mathrm{}l_n>`$, with the $`l_i`$ giving the directions of the subsequent hops, i.e. $`l_i\{\pm 1,\pm 2,\pm 3,\pm 4\}`$. So for example the terms occurring in (3) are denoted by $`<\mu ,\mu >`$ and $`<\mu ,\mu >`$, where the minus in front of the second path is due to the relative minus sign in (3).

Using this notation we can now write the most general Dirac operator on the lattice in the form:

$`D`$ $``$ $`\text{1I}[s_1+s_2{\displaystyle \underset{l_1}{}}<l_1>+s_3{\displaystyle \underset{l_2l_1}{}}<l_1,l_2>`$ (6)
$`+s_4{\displaystyle \underset{l_1}{}}<l_1,l_1>\mathrm{}]`$
$`+`$ $`{\displaystyle \underset{\mu }{}}\gamma _\mu {\displaystyle \underset{l_1=\pm \mu }{}}s(l_1)[v_1<l_1>+v_2{\displaystyle \underset{l_2\pm \mu }{}}[<l_1,l_2>`$ (8)
$`+<l_2,l_1>]+v_3<l_1,l_1>\mathrm{}]`$
$`+`$ $`{\displaystyle \underset{\mu <\nu }{}}\gamma _\mu \gamma _\nu {\displaystyle \underset{\genfrac{}{}{0pt}{}{l_1=\pm \mu }{l_2=\pm \nu }}{}}s(l_1)s(l_2){\displaystyle \underset{i,j=1}{\overset{2}{}}}ϵ_{ij}\left[t_1<l_i,l_j>\mathrm{}\right]`$ (9)
$`+`$ $`{\displaystyle \underset{\mu <\nu <\rho }{}}\gamma _\mu \gamma _\nu \gamma _\rho {\displaystyle \underset{\genfrac{}{}{0pt}{}{l_1=\pm \mu ,l_2=\pm \nu }{l_3=\pm \rho }}{}}s(l_1)s(l_2)s(l_3){\displaystyle \underset{i,j,k=1}{\overset{3}{}}}ϵ_{ijk}\times `$ (11)
$`\left[a_1<l_i,l_j,l_k>\mathrm{}\right]`$
$`+`$ $`\gamma _5{\displaystyle \underset{\genfrac{}{}{0pt}{}{l_1=\pm 1,l_2=\pm 2}{l_3=\pm 3,l_4=\pm 4}}{}}s(l_1)s(l_2)s(l_3)s(l_4){\displaystyle \underset{i,j,k,n=1}{\overset{4}{}}}ϵ_{ijkn}\times `$ (13)
$`\left[p_1<l_i,l_j,l_k,l_n>\mathrm{}\right].`$
The variables $`l_i`$ run over all of $`\{\pm 1,\pm 2,\pm 3,\pm 4\}`$ unless restricted otherwise. We use the abbreviation $`s(l_i)=`$ sign$`(l_i)`$. By $`ϵ`$ we denote the totally anti-symmetric tensors with 2,3 and 4 indices. We choose the normalization of the elements of the Clifford algebra such that the elements appear as all possible products of the $`\gamma _\mu `$ without any extra factors of $`i`$. For this normalization the symmetries C and $`\gamma _5`$-hermiticity render all coefficients $`s_i,v_i,t_i,a_i`$ and $`p_i`$ real.

The above mentioned structure of paths appearing in groups is obvious in (13) and we ordered the terms according to the length of the paths they contain. The paths in each group are related by symmetries and up to the sign factors have to come with the same real coefficient. All paths within a group have the same length. It has to be stressed, that in (13) for each element of the Clifford algebra we show only the leading terms, i.e. the terms with the shortest paths, of an infinite series of groups of paths. The dots indicate that we omitted groups with paths that are longer than the terms we show.

Let us now insert our expanded Dirac operator $`D`$ into the Ginsparg-Wilson equation. To that purpose we define:

$$E=D\gamma _5D\gamma _5+\gamma _5D\gamma _5D,$$
(14)
and finding a solution $`D`$ of the Ginsparg-Wilson equation corresponds to having $`E=0`$. Computing the linear part of $`E`$ is straightforward: When evaluating the sandwich $`\gamma _5D\gamma _5`$ we find that the terms of (13) with an odd number of $`\gamma _\mu `$, i.e. vector- and pseudovector terms, pick up a minus sign, while the other terms remain unchanged. Thus when adding the two linear terms we find that the terms with an odd number of $`\gamma _\mu `$ cancel while the other terms pick up a factor of 2.

The next step is to compute the quadratic term $`\gamma _5D\gamma _5D`$. Here we have to multiply the various terms appearing in $`D`$. Each term is made of two parts, a generator of the Clifford algebra and a group of paths. The multiplication of two of these terms proceeds in two steps: Firstly the two elements of the Clifford algebra are multiplied giving again an element of the algebra. In the second step we have to multiply the paths of our two terms. This multiplication can be noted very conveniently in our notation, where multiplication of two paths simply consists of writing the paths into one long path: $`<l_1,l_2\mathrm{}l_n>\times <l_1^{},l_2^{}\mathrm{}l_n^{}^{}>=<l_1,l_2\mathrm{}l_n,l_1^{},l_2^{}\mathrm{}l_n^{}^{}>`$. It is straightforward to establish this rule by translating back to the algebraic expression of our examples (1),(3) and performing the multiplication in this notation. It can happen that, after multiplying two paths, a hop in some direction $`l_i`$ is immediately followed by its inverse $`l_i`$. These two terms cancel each other and we find $`<l_1\mathrm{}l_{i1},l_i,l_i,l_{i+1}\mathrm{}l_n>=<l_1\mathrm{}l_{i1},l_{i+1}\mathrm{}l_n>`$. This rule is used to reduce all products of paths appearing in $`\gamma _5D\gamma _5D`$ to their true length. In a final step we decompose the sums over the product terms into groups of paths related by the symmetries in the same way as we did above when constructing the most general ansatz for $`D`$. Adding the linear and quadratic terms of (14) we end up with the following expansion for $`E`$:

$`E`$ $``$ $`\text{1I}[e_1^s+e_2^s{\displaystyle \underset{l_1}{}}<l_1>+e_3^s{\displaystyle \underset{l_2l_1}{}}<l_1,l_2>`$ (16)
$`+e_4^s{\displaystyle \underset{l_1}{}}<l_1,l_1>\mathrm{}]`$
$`+`$ $`{\displaystyle \underset{\mu }{}}\gamma _\mu {\displaystyle \underset{l_1=\pm \mu }{}}s(l_1)\left[e_1^v{\displaystyle \underset{l_2\pm \mu }{}}[<l_1,l_2><l_2,l_1>]\mathrm{}\right]`$ (17)
$`+`$ $`{\displaystyle \underset{\mu <\nu }{}}\gamma _\mu \gamma _\nu {\displaystyle \underset{\genfrac{}{}{0pt}{}{l_1=\pm \mu }{l_2=\pm \nu }}{}}s(l_1)s(l_2){\displaystyle \underset{i,j=1}{\overset{2}{}}}ϵ_{ij}\left[e_1^t<l_i,l_j>\mathrm{}\right]`$ (18)
$`+`$ $`{\displaystyle \underset{\mu <\nu <\rho }{}}\gamma _\mu \gamma _\nu \gamma _\rho {\displaystyle \underset{\genfrac{}{}{0pt}{}{l_1=\pm \mu ,l_2=\pm \nu }{l_3=\pm \rho }}{}}s(l_1)s(l_2)s(l_3){\displaystyle \underset{i,j,k=1}{\overset{3}{}}}ϵ_{ijk}\times `$ (23)
$`[{\displaystyle \underset{l_4\pm \mu ,\nu ,\rho }{}}\{e_1^a[<l_i,l_j,l_k,l_4><l_4,l_i,l_j,l_k>]`$
$`+e_2^a[<l_i,l_4,l_j,l_k><l_i,l_j,l_4,l_k>]\}`$
$`+e_3^a[<l_i,l_j,l_k,l_j><l_i,l_j,l_i,l_k>]`$
$`+e_4^a<l_i,l_j,l_k,l_i>\mathrm{}]`$
$`+`$ $`\gamma _5{\displaystyle \underset{\genfrac{}{}{0pt}{}{l_1=\pm 1,l_2=\pm 2}{l_3=\pm 3,l_4=\pm 4}}{}}s(l_1)s(l_2)s(l_3)s(l_4){\displaystyle \underset{i,j,k,n=1}{\overset{4}{}}}ϵ_{ijkn}\times `$ (25)
$`\left[e_1^p<l_i,l_j,l_k,l_n>\mathrm{}\right].`$
Again the expansion of $`E`$ is an infinite series and we display only the leading groups of paths. The coefficients $`e_i^\alpha `$ are now quadratic polynomials in the original coefficients $`s_i,v_i,t_i,a_i`$ and $`p_i`$ given by

$`e_1^s`$ $`=`$ $`2s_1+s_1^2+8s_2^2+48s_3^2+8s_4^2+8v_1^2+96v_2^2+8v_3^2`$ (27)
$`+48t_1^2+192a_1^2+384p_1^2\mathrm{},`$
$`e_2^s`$ $`=`$ $`2s_2+2s_1s_2+12s_2s_3+2s_2s_4+12v_1v_2`$ (29)
$`+2v_1v_3\mathrm{},`$
$`e_3^s`$ $`=`$ $`2s_3+2s_1s_3+s_2^2+4s_3^2+2s_3s_4+4v_2^2+2v_2v_3\mathrm{},`$ (30)
$`e_4^s`$ $`=`$ $`2s_4+2s_1s_4+s_2^2+6s_3^2v_1^26t_1^224a_1^2`$ (32)
$`48p_1^2\mathrm{},`$
$`e_1^v`$ $`=`$ $`s_2v_14s_3v_22s_4v_2s_3v_3v_3t_14v_2t_1\mathrm{},`$ (33)
$`e_1^t`$ $`=`$ $`2t_1+2s_1t_12s_4t_1v_1^24v_2^22v_2v_34t_1^2`$ (35)
$`+8v_1a_18a_1^2+16t_1p_1\mathrm{},`$
$`e_1^a`$ $`=`$ $`s_2a_1+v_2t_1v_3p_1\mathrm{},`$ (36)
$`e_2^a`$ $`=`$ $`v_2t_1\mathrm{},`$ (37)
$`e_3^a`$ $`=`$ $`s_2a_12v_2p_1\mathrm{},`$ (38)
$`e_4^a`$ $`=`$ $`2s_2a_1+2v_2t_14v_2p_1\mathrm{},`$ (39)
$`e_1^p`$ $`=`$ $`2p_1+2s_1p_12s_4p_12v_1a_1+t_1^2\mathrm{}.`$ (40)
Each coefficient $`e_i^\alpha `$ is an infinite series of terms.

The crucial step in our construction is to realize that the groups of paths together with the corresponding element of the Clifford algebra, appearing in the expansion (25), are linearly independent for arbitrary background gauge field. These groups of paths together with the products of $`\gamma _\mu `$ have to be viewed as basis elements for the expansion of $`E`$. For a solution of the Ginsparg-Wilson equation we must have $`E=0`$ and thus all coefficients $`e_i^\alpha `$ have to vanish simultaneously. Hence we have rewritten the problem of finding a solution of the Ginsparg-Wilson equation to the problem of solving the system of coupled quadratic equations (40).

Before we discuss solving (40) let us first analyze our Dirac operator $`D`$ for trivial background gauge field. In this case we can compute the Fourier transform $`\widehat{D}(p)`$ of $`D`$. For small momenta we have to implement the condition $`\widehat{D}(p)=i\mathit{}+𝒪(p^2)`$, which leads to two more equations for the coefficients, the first one makes the constant term vanish and the second one requires the slope to be equal to 1.

$`0`$ $`=`$ $`s_1+8s_2+48s_3+8s_4\mathrm{},`$ (41)
$`1`$ $`=`$ $`2v_1+24v_2+4v_3\mathrm{}.`$ (42)
These two equations have to be implemented as well, and together with (40) entirely describe Dirac operators which obey the Ginsparg-Wilson equation and have the correct behavior at small momenta.

The quadratic system of equations (40), (42) is completely equivalent to the Ginsparg-Wilson equation. The known solutions of the Ginsparg-Wilson equation, i.e. the overlap operator as well as the perfect action, correspond to two different solutions of (40) and (42). Many interesting questions, such as the existence of a connected continuum of solutions of (40), (42) can be addressed in our new framework.

In the remaining part of this letter we now explore the possibility of finding solutions of our system of quadratic equations numerically. It has to be stressed, that (40), (42) is a system of infinitely many quadratic equations for infinitely many coefficients. Thus in order to be able to apply numerical methods for solving these equations one has to introduce a cutoff, i.e. we set the coefficients for terms with longer paths to zero. This cutoff comes in a natural way for proper lattice Dirac operators, since in order to maintain universality of our lattice fermions, the coefficients $`s_i,v_i,t_i,a_i,p_i`$ have to decrease in size exponentially as the length of the corresponding paths increases.

In addition, the following qualitative argument (we use only the leading order, i.e. we ignore effects of back-folding paths) shows that the structure of the Ginsparg-Wilson equation strongly supports this exponential suppression of longer paths: When setting $`E=0`$, the Ginsparg-Wilson equation (14) turns to $`D+\gamma _5D\gamma _5=\gamma _5D\gamma _5D`$. Short paths contained in the original $`D`$ combine to longer paths when evaluating the product $`\gamma _5D\gamma _5D`$ on the right-hand side. The linear term on the left-hand side requires also these longer paths with their coefficients to be contained in $`D+\gamma _5D\gamma _5`$. The coefficients of the combined longer paths are products of the coefficients of the short paths, and since the coefficients of the shortest paths are already smaller than 1 we find an exponential suppression of the coefficients for longer paths. This behavior was confirmed numerically .

Thus the physical requirement of universality as well as the structure of the Ginsparg-Wilson equation imply the coefficients to decrease exponentially in size as the length of the paths in $`D`$ increases. We now consider a truncated expansion of $`D`$ which contains only terms with paths up to a certain length. Introducing this cutoff is equivalent to what is done when computing numerically the fixed point action using a finite parametrization for the Dirac operator. The resulting Dirac operator will only be an approximate solution of the Ginsparg-Wilson equation, but the error, i.e. the norm of $`E`$, can be made arbitrarily small by including more terms in $`D`$.

After one has introduced the cutoff discussed above, e.g. one includes only terms with a maximal length of 3, the system (40), (42) turns into a finite system. Each coefficient $`e_i^\alpha `$ in the system (40) as well as the two equations (42) now consist of only finitely many terms. We numerically solve the equations corresponding to the leading terms in the expansion of $`E`$, and so determine the coefficients $`s_i,v_i,t_i,a_i,p_i`$ in the expansion (13) of $`D`$.

This approach has been implemented and carefully tested in two dimensions . We found that adding only a few simple terms to the Dirac operator improves the chiral properties tremendously, at only a small increase of the cost of a numerical treatment. The resulting Dirac operator was tested in a dynamical simulation of the 2-flavor Schwinger model and it was found that the mass of the $`\pi `$-particle in this model can be brought down, close to the values obtained with the perfect action and the overlap operator.

We are currently implementing the quadratic system in 4 dimensions and have begun to construct approximate solutions of the Ginsparg-Wilson equation along the lines outlined in this letter. Preliminary results are very encouraging, in particular we find a considerable reduction of the additive fermion mass renormalization. These results will be presented elsewhere .

Acknowledgement: The author would like to thank Ivan Hip, Christian Lang and Uwe-Jens Wiese for active interaction during the course of this work.


text:
# Goodness-of-Fit Analysis of Radial Velocities Surveys

## 1 Introduction

Surveys of radial velocities of galaxies have played a major role in the study of the large scale structure. The analysis of such surveys has been conducted in two main directions, the mapping of the local cosmography and the estimation of the cosmological parameters (cf. Dekel 1994 for a review). The Bayesian framework provides one with very elegant and powerful tools for conducting both the mapping and parameter estimation, where the recovery of the large scale structure is done by means of the Wiener filter and the parameters are estimated by maximum likelihood (MaxLike) analysis (Zaroubi et al.  1995, hereafter ZHFL). In the case where the deviations from a homogeneous and isotropic universe constitute a Gaussian random field the Wiener filter and the MaxLike are the optimal tools for performing such an analysis (ZHFL). Indeed, the MARK III catalog of radial velocities (Willick et al.  1995, 1996, 1997a) have been recently analyzed by Wiener filtering (Zaroubi, Hoffman and Dekel 1999) and by MaxLike (Zaroubi et al.  1997). The SFI survey of da Costa et al.  (1996) has been studied by MaxLike analysis by Freudling et al.  (1999) and by Wiener filtering (Hoffman and Zaroubi, unpublished). Both surveys seem to yield similar results.

In the Bayesian MaxLike analysis one calculates the posterior probability of a model to be correct given the data (ZHFL, Vogeley and Szalay 1996). Thus the model that maximizes the likelihood function, over a given parameter (or model) space, is the most likely model in that space. The MaxLike analysis cannot guarantee, however, that the most probable model is indeed consistent with the data. It provides only a relative measure for models to be correct. It is common to adopt an independent measure for the goodness-of-fit, which is often given by the requirement that the reduced $`\chi ^2`$ is close to unity. Often, when the most likely model (given the data) passes also the goodness-of-fit test one assumes that the ‘correct’ model has been nailed down. Here, the $`\chi ^2`$ test is expanded and a much more critical test is suggested and then applied to the Mark III and SFI surveys.

The $`\chi ^2`$ ’goodness-of-fit’ is based on the assumptions that all the random variables that affect the observables are normally distributed. In the cosmological context this applies to both the underlying dynamical (e.g. density and velocity) field and the statistical errors. Thus for a survey containing $`N`$ data observables (e.g. radial velocities) the $`\chi ^2`$ of the system of $`N`$ degrees of freedom (DOF) is calculated given the model that maximizes the likelihood function. The goodness-of-fit is measured by how close is the $`\chi ^2/\mathrm{DOF}`$ to unity. This provides a global measure for the consistency of the data with the model, as it includes all the observables. A situation might occur of some ‘conspiracy’ where different parts of the data deviate from the predictions of the model, but when combined together they ‘conspire’ to yield a reasonable $`\chi ^2`$. A much stronger test on the model is to decompose the data into statistical independent eigenmodes and observe the $`\chi ^2`$ behavior of the independent modes. Eigenmode analysis, also known as principal component analysis (PCA) and the Karhunen-Loeve transform, is not a new tool in the field. It has been applied to studies of redshift surveys (Vogeley and Szalay 1996),the cosmic microwave background (Bunn 1997, Bond 1995) and more recently radial velocities surveys (Hoffman, 1999). The later study is extended here to perform the ’goodness-of-fit’ test on a mode-by-mode basis. The basic formalism is presented in § 2, and its application to the Mark III and SFI surveys is given in § 3. Our results are discussed and the conclusions are summarized in § 4.

## 2 Eigenmode Analysis of Radial Velocities

Consider a data base of radial velocities $`\{u_i\}_{i=1,\mathrm{},N}`$, where

$$u_i=𝐯(𝐫_i)\widehat{𝐫}_i+ϵ_i,$$
(1)
$`𝐯`$ is the three dimensional velocity, $`𝐫_i`$ is the position of the i-th data point and $`ϵ_i`$ is the statistical error associated with the i-th radial velocity. The assumption made here is of a cosmological model that well describes the data, that systematic errors have been properly dealt with and that the statistical errors are well understood. The data auto-covariance matrix is then written as:

$$R_{ij}<u_iu_j>=\widehat{𝐫}_j<𝐯(𝐫_i)𝐯(𝐫_j)>\widehat{𝐫}_j+\sigma {}_{ij}{}^{2}.$$
(2)
(Here $`<\mathrm{}>`$ denotes an ensemble average.) The last term is the error covariance matrix. The velocity covariance tensor that enters this equation was derived by Górski (1988, see also Zaroubi, Hoffman and Dekel 1999) and it depends on the power spectrum and cosmological parameters.

The eigenmodes of the data covariance matrix provides a natural representation of the data:

$$R\eta ^{(i)}=\lambda _i\eta ^{(i)}$$
(3)
The set of $`N`$ eigenmodes $`\{\eta ^{(i)}\}`$ constitutes an orthonormal basis and the eigenvalues $`\lambda _i`$ are arranged in decreasing order (in absolute values). A new representation of the data is given by:

$$\stackrel{~}{a}_i=\eta {}_{j}{}^{(i)}u_{j}^{}$$
(4)
This provides a statistical orthogonal representation, namely:

$$\stackrel{~}{a}_i\stackrel{~}{a}_j=\lambda _i\delta _{ij}$$
(5)
The normalized transformed variables are defined by:

$$a_i=\frac{\stackrel{~}{a}_i}{\sqrt{\lambda _i}}$$
(6)
Eq. 5 is written now as:

$$a_ia_j=\delta _{ij}$$
(7)
Note that as the modes are statistically independent one can measure the $`\chi ^2`$ of a given mode, $`\chi {}_{i}{}^{2}=a_i^2`$, and the cumulative reduced $`\chi ^2`$ is given by:,

$$\chi {}_{M}{}^{2}=\frac{1}{M}\underset{i=1}{\overset{M}{}}a_i^2$$
(8)
For normally distributed errors and a Gaussian random velocity field the $`a_i`$’s are normally distributed with zero mean and a variance of unity.

In addition the probability of finding such $`\chi _M^2`$ is calculated as well. The probability is defined by

$`P(\chi {}_{M}{}^{2})`$ $`=`$ $`P_{\chi ^2}(M\chi {}_{M}{}^{2},M)\mathrm{for}P_{\chi ^2}(M\chi {}_{M}{}^{2},M)<0.5`$ (9)
$`=`$ $`\mathrm{\hspace{0.33em}\hspace{0.33em}1}P_{\chi ^2}(M\chi {}_{M}{}^{2},M)\mathrm{otherwise},`$ (10)
where $`P_{\chi ^2}(x,M)`$ is the probability that a random variable drawn from a $`\chi ^2`$ distribution with $`M`$ degrees of freedom is less than a given value $`x`$.

## 3 Differential $`\chi ^2`$ Analysis

Here the goodness-of-fit of the Mark III and SFI surveys is studied. The models studied here are the MaxLike solutions for these surveys, which are slightly different from one another. The most likely model given Mark III is a tilted-CDM (T-CDM) of $`\mathrm{\Omega }_0=1,h=0.75`$ and $`n=0.8`$ where $`\mathrm{\Omega }_0`$ is the cosmological density parameter, $`h`$ is Hubble’s constant in units of $`100km/s/Mpc`$ and $`n`$ is the power spectrum index (Zaroubi et al.  1977). The most likely model given SFI is an open CDM (OCDM) of $`\mathrm{\Omega }_0=0.79,h=0.6`$ and $`n=0.92`$ (Fruedling et al. 1999). For both cases the MaxLike best model has a total $`\chi _{M=N}^2`$ very close to unity. Thus, from the point of view of the integral $`\chi ^2`$ the MaxLike solutions seem to be very consistent with the data. This is extended to perform a differential $`\chi ^2`$ analysis, namely to study the $`\chi ^2`$ behavior across the modes spectrum.

To study the robustness of this probe it is first applied to a linear mock catalog of Mark III, constructed from an unconstrained realization of the velocity field. This field is sampled at the location of the Mark III data points, to which normally distributed errors are added according to Mark III’s error covariance matrix. The cumulative $`\chi ^2`$ of such a catalog should oscillate around unity, given that the model used to generate the catalog is known. Indeed, this has been confirmed by an analysis of a few linear mock catalogs of Mark III. The probabilities of obtaining such $`\chi ^2`$ distribution lies comfortably within the $`90\%`$ confidence level. The non-trivial result of this test is that the very poor sampling of the long wavelength Fourier waves, i.e. cosmic variance, does not affect the goodness-of-fit test.

The differential $`\chi ^2`$and its associated probability of the Mark III and SFI surveys are presented in Fig. 1, each case analyzed in its maximum likelihood solution. A clear trend is noticed, namely over almost the entire mode spectrum the cumulative $`\chi ^2`$ increases monotonically. When all modes are included the total $`\chi ^2`$/DOF is indeed close to $`1`$, but if we had to take half the modes, starting from the top or the bottom, a very different $`\chi ^2`$ would have obtained.

The differential $`\chi ^2`$ analysis is repeated for the currently popular model of $`\mathrm{\Lambda }`$-CDM($`\mathrm{\Omega }_0=0.4,h=0.6`$ and $`n=1`$; Fig. 1). Indeed, the same trend is found in this case as well but the total $`\chi ^2`$converges to a value outside the $`90\%`$ confidence level.

The conclusions that follows is that for both data sets, Mark III and SFI, and for a variety of theoretical models the differential $`\chi ^2`$ increase monotonically with the mode number (with the exception of the first 10 modes of the Mark III). The theoretical expectation is that if indeed the data is consistent with the assumed model then $`\chi _M^2`$ will fluctuate around unity. The probability of observing such a trend given a model is very small across most of the mode number range.

## 4 Discussion

What have we learned from the differential $`\chi ^2`$ analysis? It has been found that even the most probable CDM-like model, the one that maximizes the likelihood function given the data, is not fully consistent with the data. The cumulative $`\chi ^2`$ has been calculated both downwards and upwards (namely starting from the modes with the largest and smallest eigenmode, respectively). Over more than $`90\%`$ of the modes the cumulative $`\chi ^2`$ lies well outside the $`90\%`$ confidence level, indicating a very small probability of measuring such data given the assumed model. Over most of the mode number range $`\chi _M^2`$ increases monotonically. It is this behavior of the $`\chi ^2`$ which indicates a systematic inconsistency of the model with the data. The assumed model actually contains two ingredient, the theoretical power spectrum and the error model. However, the present analysis cannot indicate which one is to be ‘blamed’ for the systematic trend. It should be noted here that apart from the first few ($`1020`$) modes there is a clear correlation of the eigenvalues with its weighted mean distance (of data points of the given mode). Namely, the variance associated with a mode (i.e. its eigenvalue) increases with its mean distance (Zehavi, private communication, Silverman et al.  in preparation). It follows that the $`\chi ^2`$ trend seen here is closely correlated with the distance and that the data ‘asks’ for less power on large scales than the model (power spectrum and noise) provides. A detailed study of the power spectrum and error model possible modifications is to be given elsewhere (Silverman et al.  in preparation). (Note that these first $`1020`$ modes are the ones dominated by the underlying velocity field and not the noise, Hoffman 1999.)

The cosmological implications of the present findings are that either the error and/or the theoretical model need to be modified. The theoretical model assumed in the analysis of large scale radial velocity surveys is that the velocities are drawn from a Gaussian random field defined by a given power spectrum. The present study might indicate the inconsistency of the power spectrum with the data. A less likely possibility is that it indicates a departure from the Gaussian statistics. Alternatively, the present work might indicate a systematic error that has not been accounted for that causes this trend. Still another possibility is that of an indication for a velocity bias.

The conclusions reached here should not be taken as a contradiction of the results of Zaroubi et al.  (1997) and Fruedling et al.  (1999), but rather as extending and complementing them. The Bayesian MaxLike analysis can be performed only within the assumed parameter/model space. The differential $`\chi ^2`$ allows one to go beyond this and analyze the nature of the agreement, or the lack of it, between a given model and the data on a mode by mode basis.

The PCA transforms the data to a statistically independent representation and enables the study of the compatibility of the data with the model on a mode by mode basis. This differential analysis is in contrast to the more traditional approach where a data set is analyzed as a whole. The differential $`\chi ^2`$ analysis should be performed together with the Bayesian MaxLike analyze and complement it. This should be useful in fields where the MaxLike is the basic tool of analysis such as the mapping of the CMB angular fluctuations and the study of redshift surveys as well as all radial velocity surveys. The present analysis can prove to be very useful and powerful in those fields where systematic errors play a crucial roles, such as redshift and radial velocities surveys.

We have benefited from many interesting discussions with Avishai Dekel, Zafrir Kolatt, Ofer Lahav, Lior Silverman, Simon White and Idit Zehavi. The hospitality of the Racah Inst. Physics and the Max Planck Institut fur Astrophysik is gratefully acknowledged. This research has been partially supported by a Binational Science Foundation grant 94-00185 and an Israel Science Foundation grant 103/98.