text:
# Effect of halo modeling on WIMP exclusion limits

## I Introduction

Arguably the best motivated non-baryonic dark matter candidate is the neutralino (the lightest supersymmetric particle), and current direct detection experiments are just reaching the sensitivity required to probe the relevant region of parameter space . The most stringent exclusion limits on Weakly Interacting Massive Particles (WIMPs) in general currently come from the Edelweiss  and Cryogenic Dark Matter Search (CDMS) experiments , with competitive constraints also having been produced by Heidelberg-Moscow  and IGEX . The sensitivity to WIMPs will be improved significantly in the short term future by the continued operation of Edelweiss, and CDMS moving in a low background environment at the Soudan mine , and in the longer term by, for instance, the planned GENIUS project .

The direct detection event rate, and its energy distribution, depend crucially on the WIMP speed distribution. Data analyzes nearly always assume a standard smooth halo model with isotropic maxwellian velocity distribution. The change in the expected signal has been calculated for various non-standard halo models of varying degrees of sophistication . For models which are effectively close to maxwellian, while there may be a significant change in the annual modulation and angular dependence of the signal, the change in the mean (averaged over time and recoil direction) differential event rate is typically small . Models with triaxiality or velocity anisotropy may however produce a significant change even in the mean differential event rate . Furthermore all of the non-standard halo models which have previously been considered are essentially smooth***An exception is Sikivie’s late infall model , which contrary to the standard picture of halo formation in CDM cosmologies (see e.g. Refs. ) assumes axial symmetry and cold collapse.. N-body simulations, however, produce dark matter halos which as well as being triaxial with anisotropic velocity distributions  also contain substructure . A number of groups have recently investigated the local dark matter distribution numerically , using different methods and reaching, to some extent, different conclusions.

Triaxiality, anisotropy and clumping in the WIMP velocity distribution could potentially have a significant effect on the WIMP direct detection signal. Constraints (and in the future possibly best fits) calculated assuming a standard maxwellian halo could be erroneous . On the other hand, more optimistically, it might be possible to derive useful information about the local velocity distribution, and hence the formation of the galactic halo, if WIMPs were detected . Belli et. al.  have recently reanalyzed the DAMA collaborations annual modulation signal  for a range of halo models, finding that the allowed region of WIMP mass–cross-section parameter space is significantly enlarged. This illustrates that it is important to take into account uncertainties in halo modeling when comparing exclusion limits and/or allowed regions from different experiments.

Given the importance of the local dark matter distribution for direct detection experiments we devote Sec. II to a detailed discussion of the global properties of real and simulated dark matter halos and recent work on the local dark matter distribution . In Sec. III we examine the effect of realistic halo modeling on exclusion limits. We first investigate triaxial and anisotropic halos models, with parameter choices motivated by the observations and simulations, and then, more speculatively examine the possible effects of small subhalos.

## II Galaxy halos

### A Global properties

Observational constraints on the structure of dark matter halos depend on the relation of luminous tracer populations to the underlying density distribution, and are complicated by galactic structure and projection effects. Ref.  concludes that in the outskirts of spiral galaxies the intermediate-to-long axis ratio is likely to be greater than $`0.8`$, while the short-to-long axis ratio is largely unconstrained with values in the range $`0.31.0`$ reported, with some correlation between the method used and the value obtained (see Refs.  and  for details and references). For the Milky Way (MW), analysis of local stellar kinematics gives an estimate for the short-to-long axis value of $`0.7\pm 0.1`$ , while the great circle tidal streams observed from the Sagittarius dwarf galaxy rule out ratios of less than 0.7 in the outer halo at high confidence  (in a flattened potential angular momentum is not conserved, so that orbits precess and tidal streams lose their coherence). The (an)isotropy of velocities in the MW halo is even harder to probe, however there is some evidence that galactic globular clusters may have preferentially radial orbits .

Given the difficulties involved in ‘observing’ galaxy halos it makes sense to turn to numerical simulations for information on their possible structure. Current simulations of galaxy halos within a cosmological context can resolve sub-kpc scales (see e.g. ). Discrepancies between the halos produced in these simulations (which have lots of surviving dwarf galaxy sized subhalos and steep central profiles) and observations, have led to claims of a crisis for the cold dark matter model (see Ref.  and references therein for extensive discussion). Most relevant for the local dark matter distribution is the subhalo problem which may be, at least partly, due to complications in comparing the observed luminous matter with the dark matter distribution from the simulations. In particular it has been argued that gas accretion onto low mass halos may be inhibited after reionization so that a large fraction of the subhalos remain dark . It has also been shown that if the observed dwarf galaxies themselves have dark halos, then their masses have been underestimated and correcting for this would go toward resolving the discrepancy with observations . The survival of subhalos is at least partly due to their concentrated profiles, so any modification to the simulations which produced halos with shallower central profiles could also reduce the number of surviving subhalos. Despite the ongoing debate regarding the detailed comparison of the small scale properties of simulated halos with observations, cosmological simulations may still provide us with useful information about the global properties of dark halos.

The shape of simulated halos varies, not just between different halos of the same mass, but also as a function of radius within a single halo, strongly if the halo has undergone a major merger relatively recently. Two high resolution Local Group halos studied in detail in Ref.  have axis ratios of $`1:0.78:0.48`$ and $`1:\mathrm{0.45.0.38}`$ at the solar radius and $`1:0.64:0.40`$ and $`1:0.87:0.67`$ at the virial radius. Adding dissipative gas to simulations tends to preserve the short-to-long axis ratio while increasing the intermediate-to-long axis ratio .

The anisotropy parameter $`\beta (r)`$, defined as

$$\beta (r)=1\frac{<v_\theta ^2>+<v_\varphi ^2>}{2<v_\mathrm{r}^2>},$$
(1)
where $`<v_\theta ^2>`$, $`<v_\varphi ^2>`$ and $`<v_\mathrm{r}^2>`$ are the means of the squares of the velocity components evaluated at radius $`r`$, also varies with radius. Typically $`\beta (r)`$ grows, although not monotonically, from roughly zero in the center of the halo to close to one at the virial radius, with non-negligible variation between halos (see Fig. 23 of Ref. ). The high resolution galactic mass halos studied in Ref.  have $`\beta (R_{})`$ in the the range 0.1-0.4, corresponding to radially biased orbits.

### B Local dark matter distribution

In CDM cosmologies structure forms hierarchically, from the top down . Small objects (often known as subhalos) form first, with larger objects being formed progressively via mergers and accretion. The internal structure of large galaxy size halos is determined by the dynamical processes which act on the accreted components. Dynamical friction causes subhalos with mass $`M10^9M_{}`$ to spiral toward the center of their parent halo within a Hubble time , while the tidal field of the main halo can strip material away from a subhalo  producing tidal streams along its orbit .

The local dark matter distribution can not be probed directly by cosmological simulations; the smallest subhalos resolvable in the highest resolution simulations have mass of order $`10^7M_{}`$ and it is not possible to fully resolve substructure within subhalos. Little substructure is found within the central regions of simulated halos, however it is not known whether the subhalos have been destroyed by tidal stripping or if this is purely a resolution effect . This is crucial for the local dark matter distribution as the solar radius ($`R_{}8`$ kpc) is small compared with the radius of the MW halo, which from observations is thought to be in excess of 100 kpc , while simulated halos with the same peak circular velocity as the MW have virial radiiThe virial radius is the radius which separates the virialized and infall regions of simulated halos and has a mean density within it of 178 (100) times the critical density at $`z=0`$ in the standard($`\mathrm{\Lambda }`$)CDM cosmology . of order 200 kpc. The computing power required to directly probe the local dark matter distribution will probably not be available for a decade or so , therefore other numerical  and semi-analytic  approaches have been used to address the problem.

Stiff, Widrow and Frieman  employ a semi-analytical approach, calculating the subhalo distribution as a function of mass and accretion redshift using the extended Press-Schechter formalism  and then evolving individual late accreting smooth subhalos within a growing spherical halo, to find the probability distribution of the over-density at the solar radius. They find that there is a high (of order one) probability that there is a density enhancement of $`3\%`$ of the mean halo density in the solar neighborhood, and the probability of an enhancement roughly equal to the background density is non-negligible (of order 0.01)These probabilities are lower limits as they only include the contributions of subhalos accreted after $`z=1`$.. They show that if the velocity of the clump, or stream, with respect to the Earth is high enough it will produce a shoulder in the differential event rate at high energies.

Moore et. al.  take a region with dynamical properties similar to the Local Group, resimulated at higher resolution from a standard CDM cosmological simulation, and identify a subhalo with mass similar to the Draco dwarf galaxy. They then resimulate this subhalo up until its merger with the parent halo, at which point roughly 10$`\%`$ of its mass is in the form of (sub)subhalos. To assess the effect of tidal friction, subhalos with smooth density profiles are evolved within a smooth Galactic potential and it is found that small subhalos orbiting at the solar radius which are accreted onto the Galactic halo early, retain most of their mass due to their high central densities.

Moore et. al. conclude that the phase space distribution at the solar radius will depend crucially on the Galaxy’s merger history and on the internal structure of the smallest subhalos, arguing that it is possible that the local dark matter density could be zero or that a single dark matter stream with small velocity dispersion could dominate or that many tidal streams could overlap to give a smooth distribution. The solution to this problem depends on the extent to which the substructure within subhalos is destroyed prior to their accretion onto the main halo. The subhalos are of course also formed hierarchically themselves from smaller subhalos. The free-streaming length of neutalinos is so small that the first clumps to form have mass $`M10^{12}M_{}`$, however, and to follow the accretion and destruction of subhalos though such a large hierarchy of scales would be decidedly non-trivial.

Helmi, White and Springel  take a cosmological $`\mathrm{\Lambda }`$CDM simulation, where the second largest cluster has been resimulated at higher resolution and then scaled down in size to match the MW. Most of the mass in the inner halo has been in place for 10 Gyr and the smallest resolvable subhalos do not survive to the present day at the solar radius. They argue that the only way a small subhalo with high density could exist in the solar neighborhood at $`z=0`$ is if it were first accreted by a large halo which subsequently had a major merger with the main progenitor of the Galactic halo. They conclude that the local dark matter velocity distribution is well approximated by a smooth multi-variate gaussian, with clumps of high velocity particles present if the MW halo has undergone a recent major merger.

In summary numerical simulations produce galaxy halos which are significantly triaxial and anisotropic, with the shape and anisotropy of a halo depending on its individual merger history, a picture which is broadly supported by observations. This indicates that, even if the local velocity distribution is relatively smooth, the standard spherical isotropic maxwellian halo model may not be a good approximation. Furthermore galaxy halos are formed hierarchically from the accretion of smaller subhalos, which may not be completely destroyed by tidal friction, so that the local dark matter could be distinctly non-smooth. It is even possible that the dark matter could be distributed largely in small dense clumps. There is currently no consensus on the local dark matter velocity distribution however, with the results obtained depending on the method used to extrapolate to small scales below the resolution limit of cosmological simulations, so this possibility should be regarded as speculative.

## III Effects on experimental analysis

The differential WIMP event rate due to scalar interactions can be written in terms of the WIMP scattering cross section on the proton, $`\sigma _\mathrm{p}`$ :

$$\frac{\mathrm{d}R}{\mathrm{d}E}=\zeta \sigma _\mathrm{p}\left[\frac{\rho _{0.3}}{\sqrt{\pi }v_0}\frac{(m_\mathrm{p}+m_\chi )^2}{m_\mathrm{p}^2m_\chi ^3}A^2T(E)F^2(q)\right],$$
(2)
where the local WIMP density, $`\rho _\chi `$ is normalized to a fiducial value $`\rho _{0.3}=0.3\mathrm{GeV}/\mathrm{cm}^3`$, such that $`\zeta =\rho _\chi /\rho _{0.3}`$, $`m_A`$ is the atomic mass of the target nuclei, $`E`$ is the recoil energy of the detector nucleus, and $`T(E)`$ is defined as

$$T(E)=\frac{\sqrt{\pi }v_0}{2}_{v_{\mathrm{min}}}^{\mathrm{}}\frac{f_v}{v}dv,$$
(3)
where $`f_v`$ is the WIMP speed distribution in the rest frame of the detector, normalized to unity, and $`v_{\mathrm{min}}`$ is the minimum detectable WIMP speed

$$v_{\mathrm{min}}=\left(\frac{E(m_\chi +m_A)^2}{2m_\chi ^2m_A}\right)^{1/2}.$$
(4)
### A Triaxiality and anisotropy

We will first examine the change in the WIMP speed distribution for triaxial and anisotropic, but still smooth, halo models. To date two self-consistent triaxial and/or anisotropic halo models have been studied in relation to WIMP direct detection: the logarithmic ellipsoidal model  and the Osipkov-Merritt anisotropy model , studied in Ref. . We will extend the previous work by focusing on parameters which span the range of halo properties discussed in Sec. II above and including the detector response.

The logarithmic ellipsoidal model  is the simplest triaxial generalization of the isothermal sphere and the velocity distribution can be approximated by a multi-variate gaussian on either the long or the intermediate axis (see Appendix A for further details)§§§Of course there is no reason to expect the Sun to be located on one of the axes of the halo.. We consider parameter values $`p=0.9,q=0.8`$ corresponding to axial ratios $`1:0.78:0.48`$ (as used in Ref. ), and $`p=0.72,q=0.70`$ corresponding to $`1:0.45:0.38`$, and locations on the long and intermediate axes. The first set of axial ratios is typical of the values found in simulations and roughly consistent with observations, while the second is arguably rather extreme. We use values of the anisotropy parameter $`\gamma `$ (which in the spherical limit $`p=q=1`$ is related to $`\beta `$ by $`\gamma =2\beta `$) which give $`\beta =0.1`$ and $`0.4`$, and are tabulated in Table I. The speed distributions, in the rest frame of the Sun normalized to unity, are plotted in Fig. 1 along with that for the standard maxwellian halo model. For both positions the triaxial models have a wider spread in speeds than the standard model, so that the differential event rate will decrease less rapidly with increasing recoil energy, but the change is small on the major axis. This is because the change in the speed distribution is largely determined by the velocity dispersion in the $`\varphi `$ direction. On the major axis, for parameter values which give $`0.1<\beta <0.4`$, all three components of the velocity have roughly the same dispersion, whereas on the intermediate axis the velocity dispersion in the $`\varphi `$ direction is significantly larger than that in the $`z`$ direction. Note that the speed distributions we consider deviate less from the standard maxwellian than those considered by Refs.  and  as they use more extreme values for $`\gamma `$ (namely 16 and -1.78).

In the Osipkov-Merritt (OM) model (see Appendix B for further details), which assumes a spherically symmetric density profile, the velocity anisotropy varies as a function of radius as

$$\beta (r)=\frac{r^2}{r^2+r_\mathrm{a}^2},$$
(5)
so that the degree of anisotropy increases with increasing radius, as is found in numerical simulations. Following Ref.  we assume a NFW  density profileVarying the inner slope of the density profile does not significantly affect the local velocity distribution . with scale radius $`r_\mathrm{s}=20`$ kpc. We use values of the anisotropy radius $`r_\mathrm{a}=20,12,`$ and $`9.8`$ kpc which correspond to $`\beta (R_0)=0.14,0.31`$ and $`0.4`$ respectively. For the first two values analytic fitting functions for the distribution function have been provided by Widrow . The resulting speed distributions are plotted in Fig. 2 along with that for the standard maxwellian halo model. The excess at large $`v`$ is due to the increased number of particles on very elongated, nearly radial orbits .

To assess the effect of changes in the speed distribution on exclusion limits we need to take into account the detector response, including the difference between the observed energy of an event and the actual recoil energy, non-zero energy threshold and energy resolution (see Ref.  for further details), as these factors may blur out the effects of changes in the speed distribution. We consider a $`\mathrm{Ge}^{76}`$ detector with energy threshold $`E_\mathrm{T}=4`$ keV with the same properties (resolution and form factor) as that used by the IGEX experiment , which is optimized for detecting double-beta decay. The resulting differential event rates, per kg per day per keV, are plotted in Fig. 3, for the OM speed distributions in Fig. 2, for WIMPs with mass $`m_\chi =50`$ GeV and cross-section $`\sigma _\mathrm{p}=10^{45}\mathrm{m}^2`$. We see that, in this model, the differential event rate does not deviate linearly from that of the standard halo model as the degree of anisotropy is increased. Since future detectors, optimized for WIMP detection, will have lower thresholds and better energy resolution, we also plot the differential event rates for an ‘ideal’ (i.e. completely unrealistic) $`A=76`$ detector where the full recoil energy is detected, the energy resolution is perfect and $`E_\mathrm{T}=0`$ keV Somewhat counter-intuitively, at low energies the differential event rate for the IGEX detector is higher than that for the ideal detector, due to the finite energy resolution of the IGEX detector.. The difference between the differential event rates is then largest at small recoil energies, and would therefore be more significant for an experiment with a lower threshold energy.

In Fig. 4 we plot the exclusion limits found from the IGEX data by requiring that the data in no more than one energy bin exceeds its 99.77$`\%`$ confidence limit, so as to produce 90$`\%`$ overall confidence limits , for the logarithmic ellipsoidal model and for the OM anisotropy model. We also plot the exclusion limits from the Heidelberg-Moscow (HM) experiment  for the OM anisotropy model in Fig. 5. Comparing Figs. 4 and 5 we see that the change in the exclusion limits depends not only on the halo model under consideration, but also on the data being used; for IGEX the change in the exclusion limits is largest for large WIMP masses, while for HM the change is largest for small WIMP masses. For different WIMP masses, different energy ranges can be most constraining; for the IGEX data the lowest energy bin is always the most constraining, while for HM as the WIMP mass increases the constraint comes from higher energy bins. It should therefore be borne in mind when comparing exclusion limits from different experiments, that changing the assumed WIMP speed distribution will affect the limits from different experiments differently.

The change in the exclusion limits is not huge (of order tens of per-cent) for the experiments we have considered, however these experiments are not optimized for WIMP detection. As illustrated in Fig. 3 the change in the differential event rate, and hence the exclusion limit, would be significantly larger for an experiment with better energy resolution and lower threshold energy. We have also seen that different models with the same value for the anisotropy parameter $`\beta `$ have very different speed distributions, and hence a different effect on the exclusion limits. Furthermore it is conceivable that the local WIMP velocity distribution may deviate even further from the standard maxwellian distribution, than the models that we have considered.

### B Clumps

Even if dynamical processes produce a smooth background dark matter distribution, late accreting clumps may lead to a local density enhancement and velocity clumping , and produce a shoulder in the differential event rate, if their density and velocity with respect to the earth are large enough. For the experiments we have been considering the lower energy bins are most constraining, so that only very rare high density and velocity streams would have a non-negligible effect on the exclusion limits. The effect of these late accreting clumps on the annual modulation and directional signals would be more significant however .

We will now turn our attention to the consequences of the more speculative possibility that small subhalos may survive at the solar radius. We could then be located within a subhalo with local density in excess of the mean value of $`0.3\mathrm{GeVcm}^3`$, on the other hand it is even possible that we could be in a region between clumps and streams where the WIMP density is zero . In the latter case all attempts at WIMP direct detection would be doomed to failure, and exclusion limits would tell us nothing about the WIMP cross-section. At the other extreme a tiny subhalo at the earth’s location would produce a distinctive signal and, due to the enhanced density, make it easier to detect WIMPs of a given cross-section. Subhalos with $`M10^9M_{}`$ would have negligible velocity dispersion and hence a delta-function speed distribution. The resulting theoretical differential event rate would be a step function with amplitude inversely proportional to the speed of the subhalo with respect to the earth, and position increasing with increasing relative speed and WIMP mass. Consequently for small subhalo velocities and WIMP masses there would be no constraint on the WIMP cross-section (no WIMPs would have sufficient energy to cause an observable recoil), but as the WIMP mass is increased the constraints would become much tighter as then all the WIMPs would be energetic enough to cause events of a given recoil energy. In Fig. 6 we plot the exclusion limits on $`\zeta \sigma _\mathrm{p}`$ from the IGEX data for subhalos with various relative speeds\**\**\**In reality the subhalo speed would be another unknown variable. dominating the local WIMP distribution. Unlike the smooth halo models for this data set, as the WIMP mass is increased higher energy bins are most constraining, and this leads to sharp changes in the exclusion limits due to the sharp transition from no WIMPs to all WIMPs having sufficient to energy cause events of a given recoil energy. Note that the high density of a subhalo would lead to $`\xi >1`$, so that the exclusion limits on the WIMP cross-section would be tighter (by some unknown factor) than for the standard halo model.

How dense could small scale clumps be? Dark matter clump densities are usually parameterized in terms of their concentration, $`c`$, defined as

$$c=\frac{R_{\mathrm{vir}}}{r_\mathrm{s}},$$
(6)
where $`R_{\mathrm{vir}}`$ is the virial radius and $`r_\mathrm{s}`$ is the scale radius, the radius at which the effective logarithmic slope of the density profile is equal to -2 (see e.g. Refs. ). Low mass halos typically form earlier, when the density of the Universe is higher, so that their concentrations are higher than those of larger halos  and Bullock et. al.  and Eke, Navarro and Steinmetz (ENS)  have recently constructed toy models which reproduce the scaling of concentration with mass and redshift for the simulated galaxy-sized halos.

Since the first neutralino clumps to form have mass of order $`10^{12}M_{}`$ , the first generation of halos which are formed in the same way as galaxy sized halos (via the accretion and merger of smaller clumps) will have mass of order $`10^{10}M_{}`$, and virial radii of order $`0.01`$ pc. While these clumps are tiny by cosmological standards, they are still relatively large compared to the distance traveled by the Sun in a year ($`10^4`$ pc). Extrapolating the Bullock et. al. and ENS toy models, way beyond their intended range of applicability, produces significantly different ‘guesstimates’ for $`c`$ for $`M10^{10}M_{}`$ halos of 150 and 30 respectively . Subhalos within larger halos may in fact be more concentrated than isolated halos of the same mass as dense regions tend to collapse earlier and tidal stripping may steepen their density profiles . The tidal radius  of a $`10^{10}M_{}`$ subhalo orbiting in the MW halo at the solar radius is of order $`10^3`$ pc, so the dense central regions of these subhalos should survive tidal stripping . The density of simulated halos diverges toward their center, so as a measure of the typical density of a subhalo we use the mean density within the scale radius, assuming a NFW density profile :

$$\rho (r)\frac{r_\mathrm{s}}{r(1+r/r_\mathrm{s})^2}.$$
(7)
For $`c=30(150)`$ this gives $`\rho =0.9(70)\mathrm{GeVcm}^3`$, or equivalently $`\zeta =3(200)`$. We emphasize that this rough calculation relies on the extrapolation of the concentration-mass scaling models far beyond the regime in which they have been numerically tested (isolated galaxy mass halos) and the numbers obtained should probably not be taken seriously. It does illustrate, however, that if small subhalos survive then they may lead to a significant local enhancement in the WIMP density, and that understanding the small scale structure of galactic halos is crucial for realistic modeling of the WIMP velocity distribution.

## IV Discussion

Rapid progress is being made in the field of WIMP direct detection, with experiments closing in on the sensitivity required to detect neutralinos, if they constitute a non-negligible fraction (greater than $`10^4`$) of the halo density . Data analyzes usually assume the simplest halo model: an isothermal sphere with maxwellian velocity distribution. There is no clear justification, either observational or theoretical, for this assumption apart from simplicity. In fact numerical simulations  and observations  suggest that galaxy halos are triaxial and anisotropic. The local density distribution may also be non-smooth, with late accreting subhalos leading to velocity clumping . More speculatively it is even possible that the dark matter could be mainly distributed in tiny dense clumps, so that the local density distribution could be dominated by a single clump, or could even be zero . It is therefore crucial to examine the effect of realistic halo modeling on the WIMP direct detection signal.

In this paper we have investigated the change in exclusion limits due to triaxiality, velocity anisotropy and small scale clumping, taking into account detector performance and using parameter values motivated by numerical simulations and observations. Triaxiality and velocity anisotropy lead to non-negligible changes in the exclusion limits, even for detectors with relatively poor energy resolution. Furthermore the changes are different for different data sets and depend on how the anisotropy is modeled. If the local WIMP distribution is dominated by small scale clumps then the local density may be zero (making it impossible to detect WIMPs) or significantly enhanced (making it easier to detect WIMPs with a given cross-section), and the exclusion limits are changed dramatically. Clearly the survival of subhalos at the solar radius is a very important issue for WIMP direct detection.

Even if the local WIMP distribution is smooth, to derive reliable constraints on WIMP parameters and compare results for different experiments a framework needs to be developed for dealing with the uncertainty in the WIMP velocity distribution; either the development of a framework for parameterizing deviations from a baseline model, or the establishment of an agreed set of benchmark models, spanning the the range of plausible WIMP velocity distributions.

## Acknowledgments

A.M.G. was supported by the Swedish Research Council and would like to thank Amina Helmi and Ben Moore for answering questions about their work.

## A Logarithmic ellipsoidal model

The logarithmic ellipsoidal model  is the simplest triaxial generalization of the isothermal sphere and the velocity distribution can be approximated by a multi-variate gaussian on either the long or the intermediate axis:

$$f(v)=\frac{1}{(2\pi )^{3/2}\sigma _\mathrm{r}\sigma _\varphi \sigma _\mathrm{z}}\mathrm{exp}\left(\frac{v_\mathrm{r}^2}{\sigma _\mathrm{r}^2}\frac{v_\varphi ^2}{\sigma _\varphi ^2}\frac{v_\mathrm{z}^2}{\sigma _\mathrm{z}^2}\right).$$
(A1)
On the intermediate axis

$`\sigma _\mathrm{r}^2`$ $`=`$ $`{\displaystyle \frac{v_0^2p^4}{(2+\gamma )(1p^2+q^2)}},`$ (A2)
$`\sigma _\varphi ^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2q^2p^2)}{2(1p^2+q^2)}},`$ (A3)
$`\sigma _\mathrm{z}^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2p^2)}{2(1p^2+q^2)}},`$ (A4)
and on the major axis

$`\sigma _\mathrm{r}^2`$ $`=`$ $`{\displaystyle \frac{v_0^2}{(2+\gamma )(p^2+q^21)}},`$ (A5)
$`\sigma _\varphi ^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2q^21)}{2(p^2+q^21)}},`$ (A6)
$`\sigma _\mathrm{z}^2`$ $`=`$ $`{\displaystyle \frac{v_0^2(2p^21)}{2(p^2+q^21)}},`$ (A7)
where $`p`$ and $`q`$ are constants which satisfy $`q^2p^21`$ and are related to the axial ratios of the density distribution, $`I_{1,2}`$, by

$`I_1^2`$ $`=`$ $`{\displaystyle \frac{p^2(p^2q^2+p^2q^2)}{q^2+p^2p^2q^2}},`$ (A8)
$`I_2^2`$ $`=`$ $`{\displaystyle \frac{q^2(p^2q^2p^2+q^2)}{q^2+p^2p^2q^2}}.`$ (A9)
and $`\gamma `$ is a constant isotropy parameter, which in the spherical limit $`p=q=1`$ is related to $`\beta `$ (as defined in eq. (1)) by $`\gamma =2\beta `$.

## B Osipkov-Merritt model

The distribution function of a self-gravitating system with spherically symmetric density profile $`\rho (r)`$ is given, assuming an isotropic velocity distribution, by Eddington’s formula  (see also Refs. ):

$$f(\epsilon )=\frac{1}{\sqrt{8}\pi ^2}\left[_0^\epsilon \frac{\mathrm{d}^2\rho }{\mathrm{d}\mathrm{\Psi }^2}\frac{\mathrm{d}\mathrm{\Psi }}{\sqrt{\epsilon \mathrm{\Psi }}}+\frac{1}{\sqrt{\epsilon }}\left(\frac{\mathrm{d}\rho }{\mathrm{d}\mathrm{\Psi }}\right)_{\mathrm{\Psi }=0}\right],$$
(B1)
where $`\mathrm{\Psi }(r)=\mathrm{\Phi }(r)+\mathrm{\Phi }(r=\mathrm{})`$, $`\mathrm{\Phi }(r)`$ is the potential of the system, $`\epsilon =E+\mathrm{\Phi }(r=\mathrm{})=E_{\mathrm{kin}}+\mathrm{\Psi }(r)`$, and $`E`$ and $`E_{\mathrm{kin}}`$ are the total and kinetic energy respectively. In the Osipkov-Merritt model  the distribution function also depends on the angular momentum, $`L`$, of the system through the variable $`𝒬`$:

$$𝒬\epsilon \frac{L^2}{2r_\mathrm{a}^2},$$
(B2)
where $`r_\mathrm{a}`$ is the anisotropy radius, which is related to $`\beta `$ by

$$\beta (r)=\frac{r^2}{r^2+r_\mathrm{a}^2},$$
(B3)
so that the degree of anisotropy increases with increasing radius. The distribution function in the Osipkov-Merritt model is found by replacing $`\epsilon `$ by $`𝒬`$ and $`\rho (r)`$ by $`(1+r^2/r_\mathrm{a}^2)\rho (r)`$ in eq. (B1), which then has to be solved numerically. Analytic fitting functions for have been provided by Widrow  for selected values of $`r_\mathrm{a}`$, for the NFW density profile . Note that physical models only exist for $`r_\mathrm{a}>r_{\mathrm{a},\mathrm{min}}`$, where $`r_{\mathrm{a},\mathrm{min}}`$ depends on the potential.


format: text

text:
# Cosmological Constraints from Compact Radio Source Angular Size versus Redshift Data

## 1 Introduction

Cosmological models now under consideration have a number of adjustable parameters. A simple way to determine whether a model provides a useful approximation to reality is to use many different cosmological tests to set constraints on cosmological-model-parameter values and to check if these constraints are mutually consistent (see, e.g., Maor et al. 2002; Wasserman 2002).

During the past few years much attention has been focussed on the Type Ia supernova apparent magnitude versus redshift test (see, e.g., Riess et al. 1998; Perlmutter et al. 1999; Podariu & Ratra 2000; Waga & Frieman 2000; Gott et al. 2001; Leibundgut 2001).111 The proposed SNAP satellite should provide much tighter constraints on cosmological parameters from this test (see, e.g., Podariu, Nugent, & Ratra 2001a; Weller & Albrecht 2002; Wang & Lovelace 2001; Gerke & Efstathiou 2002; Eriksson & Amanullah 2002). This cosmological test indicates that the energy density of the current universe is dominated by a cosmological constant, $`\mathrm{\Lambda }`$, or a term in the material stress-energy tensor that only varies slowly with time and space and so behaves like $`\mathrm{\Lambda }`$.

In conjunction with dynamical estimates which indicate a low non-relativistic matter density parameter $`\mathrm{\Omega }_0`$ (see, e.g., Peebles 1993), cosmic microwave background anisotropy measurements also suggest the presence of $`\mathrm{\Lambda }`$ or a $`\mathrm{\Lambda }`$-like term (see, e.g., Podariu et al. 2001b; Wang, Tegmark, & Zaldarriaga 2002; Baccigalupi et al. 2002; Durrer, Novosyadlyj, & Apunevych 2001; Scott et al. 2002; Mason et al. 2002).

However, the observed rate of multiple images of radio sources or quasars, produced by gravitational lensing by foreground galaxies, appears to favor a smaller value for $`\mathrm{\Lambda }`$ (see, e.g., Ratra & Quillen 1992; Helbig et al. 1999; Waga & Frieman 2000; Ng & Wiltshire 2001) than is indicated by the observations mentioned above. It is therefore of interest to examine the entrails of other cosmological tests.

In this paper we consider the redshift-angular size test, using the Gurvits, Kellermann, & Frey (1999) compact radio source measurements. The redshift-angular size relation is measured, for structures a few orders of magnitude larger than those considered by Gurvits et al. (1999), by Buchalter et al. (1998) for quasars, and by Guerra, Daly, & Wan (2000) for radio galaxies; we do not use these data sets in our analysis here. Vishwakarma (2001), Lima & Alcaniz (2002), and references therein, use the Gurvits et al.(1999) data to set constraints on cosmological parameters; our results are consistent with, but as discussed next extend, their analyzes.222 The Gurvits et al. (1999) data augments that of Kellermann (1993). Stelmach (1994), Stepanas & Saha (1995), Jackson & Dodgson (1996), and Kayser (1995) discuss the Kellermann (1993) data.

Cosmological applications of the redshift-angular size test require knowledge of the linear size of the “standard rod” used. Some earlier analyzes of the Gurvits et al. (1999) data appear to assume that this linear size will be determined using additional data, and so quote limits on cosmological parameters (such as $`\mathrm{\Omega }_0`$ or the cosmological constant density parameter $`\mathrm{\Omega }_\mathrm{\Lambda }`$) for a range of values of this linear size. Here we note that it is best to treat this linear size as a “nuisance” parameter (for the cosmologically relevant part of this test), that is also determined by the redshift-angular size data, and so marginalize over it (using a prior to incorporate other information about it, if needed).333 The situation here is similar to that for the redshift-magnitude test (e.g., Riess et al. 1998; Perlmutter et al. 1999) where one must marginalize over the magnitude of the standard candle used, treating it as a nuisance parameter. In fact, Gurvits et al. (1999) determine the linear size from the redshift-angular size data by using the model of Gurvits (1994).

In $`\mathrm{\S }`$ 2 we summarize our computation. Results are presented and discussed in $`\mathrm{\S }`$ 3. We conclude in $`\mathrm{\S }`$ 4.

## 2 Computation

For our analyzes here we use the redshift-angular size data of Fig. 10 of Gurvits et al. (1999), which are binned redshift-angular size data derived from measurements of 145 sources. These measurements are combined in twelve redshift bins, with about the same number of sources per bin, with the lowest and highest redshift bins centered at redshifts $`z=0.52`$ and $`z=3.6`$.

We consider two cosmological models as well as a currently popular parametrization of dark energy. These are low-density cold dark matter (CDM) dominated cases, consistent with current observational indications. The first model is parametrized by two “cosmological” parameters, $`\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ (in addition to all the other usual parameters). This model includes, as special cases, two “one parameter” models: the currently popular $`\mathrm{\Lambda }`$CDM case with flat spatial hypersurfaces and $`\mathrm{\Lambda }>0`$ (see, e.g., Peebles 1984; Efstathiou, Sutherland, & Maddox 1990; Stompor, Górski, & Banday 1995; Ratra et al. 1997; Sahni & Starobinsky 2000) and a model with open spatial hypersurfaces and no $`\mathrm{\Lambda }`$ (see, e.g., Gott 1982; Ratra & Peebles 1995; Cole et al. 1997).

We also derive constraints on the parameters of a spatially-flat model with a dark energy scalar field ($`\varphi `$) with scalar field potential energy density $`V(\varphi )`$ that at low $`z`$ is $`\varphi ^\alpha `$, $`\alpha >0`$. The energy density of the scalar field decreases with time, behaving like a time-variable $`\mathrm{\Lambda }`$ (see, e.g., Peebles & Ratra 1988; Ratra & Peebles 1988; Steinhardt 1999; Brax, Martin, & Riazuelo 2000; Carroll 2001).

In linear perturbation theory, a scalar field is mathematically equivalent to a fluid with time-dependent equation of state parameter $`w=p/\rho `$ and speed of sound squared $`c_s^2=dp/d\rho `$, where $`p`$ is the pressure and $`\rho `$ the energy density (see, e.g., Ratra 1991). The XCDM parametrization for dark energy approximates $`w`$ as a constant (see, e.g., Steinhardt 1999; Sahni & Starobinsky 2000; Carroll 2001; Huterer & Turner 2001), which is accurate in the radiation and matter dominated epochs but not in the current, dark energy scalar field dominated epoch. Nevertheless the XCDM parametrization is recommended by its simplicity, so we also determine redshift-angular size constraints on its parameters.

We want to determine how well the Gurvits et al. (1999) redshift-angular size measurements distinguish between different cosmological-model-parameter values. To do this we pick one of the above models or the XCDM parametrization, and a range of model-parameter values and compute the angular size distance $`r(z)`$ (the coordinate position for the object considered at redshift $`z`$ with the observer at the origin) for a grid of model-parameter values that span this range. An object of physical or proper length $`l`$ transverse to the line of sight subtends an angle $`\theta (z)=l(1+z)/[a_0r(z)]`$, where $`a_0`$ is the current value of the scale factor (Peebles 1993, $`\mathrm{\S }`$ 13). To determine the probability distribution of the cosmological model parameters ($`P`$), we compute

$$\chi ^2(l,P)=\underset{i=1}{\overset{12}{}}\left[\frac{\theta (l,P,z_i)\theta _{\mathrm{obs}}(z_i)}{\sigma (z_i)}\right]^2,$$
(1)
where $`\theta _{\mathrm{obs}}(z_i)`$ and $`\sigma (z_i)`$ are the observed angles and errors for each of the twelve redshift bins centered at redshifts $`z_i`$ of the Gurvits et al. (1999) data. $`P`$ represents the model parameters, for instance $`\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ in the general two-dimensional constant $`\mathrm{\Lambda }`$ case. This representation (eq. ) is exact for the case where the correlated errors between redshift bins are negligible.

With a uniform prior for the physical length $`l`$, the probability distribution (likelihood) of the cosmological model parameters is

$$L(P)=𝑑le^{\chi ^2(l,P)/2},$$
(2)
where the integral is over a large enough range of $`l`$ to include almost all of the probability. We typically integrate $`l`$ over the range 1 $`h^1`$ to 60 $`h^1`$ pc,444 Here the Hubble constant $`H_0=100h`$ km s-1 Mpc-1. sampled at 60 points, although the range depends somewhat on the model under analysis. For the constant $`\mathrm{\Lambda }`$ model we compute over the ranges $`0\mathrm{\Omega }_01`$ and $`1\mathrm{\Omega }_\mathrm{\Lambda }1`$, both sampled at 201 points; for the scalar field dark energy model we compute over the ranges $`0.005\mathrm{\Omega }_00.995`$ and $`0\alpha 8`$ sampled at 199 and 161 points, respectively; while for the XCDM parametrization we compute over the ranges $`0\mathrm{\Omega }_01`$ and $`1w0`$, both sampled at 201 points. The confidence limits are computed from the distribution of eq. (2), and we typically show 1, 2, and 3 $`\sigma `$ confidence contours which correspond to enclosed probabilities of 68.27, 95.45, and 99.73 $`\%`$, respectively.

Since $`\mathrm{\Omega }_0`$ is a positive quantity, we also consider the non-informative or logarithmic prior $`p(\mathrm{\Omega }_0)1/\mathrm{\Omega }_0`$ (Berger 1985; Gott et al. 2001) and compute confidence regions for the probability distribution $`L(P)/\mathrm{\Omega }_0`$, where $`L(P)`$ is given in eq. (2).

## 3 Results and Discussion

Figure 1 shows the Gurvits et al. (1999) redshift-angular size constraints on the general two-dimensional constant $`\mathrm{\Lambda }`$ case. Apparently these have not previously been published. These constraints are consistent with, but mostly not as constraining as those from Type Ia supernova redshift-magnitude data, except at larger values of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`\mathrm{\Omega }_0`$ (see, e.g., Podariu & Ratra 2000, Fig. 5).

Figure 2 shows the redshift-angular size constraints on the XCDM parameters. Lima & Alcaniz (2002, Fig. 2) show related constraints computed at fixed physical length $`l`$ for the same Gurvits et al. (1999) data. While the shapes are similar, the Lima & Alcaniz (2002) contours are much more constraining than those found here, largely because our procedure of marginalizing over the physical length $`l`$ also accounts for the uncertainty in the determination of $`l`$. Podariu & Ratra (2000, Fig. 2) show corresponding constraints on the XCDM parameters from the Type Ia supernova redshift-magnitude data, which are significantly more constraining that those shown in Fig. 2 here.

Figure 3 shows the constraints on the dark energy scalar field model with potential energy density $`V(\varphi )\varphi ^\alpha `$, $`\alpha >0`$ (Peebles & Ratra 1988). They are consistent with, but not as constraining as those from the Type Ia supernova redshift-magnitude data (Podariu & Ratra 2000; Waga & Frieman 2000).

## 4 Conclusion

Constraints on cosmological model parameters derived from the redshift-angular size compact radio source data of Gurvits et al. (1999) are consistent with but less constraining than those derived from the redshift-magnitude Type Ia supernova data of Riess et al. (1998) and Perlmutter et al. (1999).

Higher quality redshift-angular size data will more significantly constrain cosmological models, and in combination with high quality redshift-magnitude data will provide a check of conventional general relativity on cosmological length scales.

We are indebted to L. Gurvits for providing the binned redshift-angular size data. We acknowledge helpful discussions with J. Alcaniz, R. Daly, J. Lima, and J. Peebles, and support from NSF CAREER grant AST-9875031.


text:
# Beyond the quark model of hadrons from lattice QCD

## 1 Introduction

Lattice QCD is a first-principles approach to solving QCD non-perturbatively. It enables theorists to explore the consequences of QCD for quarks of different mass and so complements experimental studies. It is especially valuable for validitating phenomenological models. Direct comparison with experiment is hampered by two obstructions: (i) the continuum limit (lattice spacing $`a0`$) should be taken but this is not computationally feasible in many cases and (ii) the $`u`$ and $`d`$ quark masses are too light to be successfully explored on a lattice so that one has to rely on chiral perturbation theory and extrapolate from heavier quark masses. Nevertheless, lattice QCD enables exploration of the hadron spectrum and matrix elements and can be used to quantify the success (and failure) of quark models. One area where the quark model is deficient is in flavour singlet mesons. Here the gluonic degrees of freedom are likely to make a significant contribution.

## 2 Flavour singlet mesons

In lattice QCD, one studies mesons by creating a quark anti-quark pair at time 0 with the quantum numbers of that meson and then annihilating the quark anti-quark pair at time $`t`$. The quark sources and sinks are then combined using the fully non-perturbative quark propagators determined on the lattice. For flavour singlet mesons, there will be two contributions: disconnected $`D(t)`$ and connected $`C(t)`$ as illustrated in fig. 1.

At large t where ground state contributions dominate these measured correlations satisfy $`C(t)=ce^{m_1t}`$ and $`C(t)+D(t)=de^{m_0t}`$ where $`m_0`$ is the flavour singlet mass and $`m_1`$ the flavour non-singlet mass. Then by a study of $`D/C`$ which is given by $`(d/c)\mathrm{exp}((m_1m_0)t)1`$ one can explore the mass splitting between flavour singlet and non-singlet which is a measure of the OZI violating gluonic contribution arising from the disconnected diagram.

From a lattice study of this for all mesonic quantum numbers all-ozi  the only significant disconnected contributions were found for scalar and for pseudoscalar mesons. This is not unexpected: both from the phenomenology of the observed meson masses and from the expectation that the lightest glueball has scalar quantum numbers and that topological charge fluctuations have pseudoscalar quantum numbers.

For a review of the situation concerning pseudoscalar mesons see ref eta . We now discuss the glueballs and their impact on the meson spectrum.

## 3 Glueballs and scalar mesons

Glueballs are defined to be hadronic states made primarily from gluons. The full non-perturbative gluonic interaction is included in quenched QCD. In the quenched approximation, there is no mixing between such glueballs and quark - antiquark mesons. A study of the glueball spectrum in quenched QCD is thus of great theoretical value.

This has been studied extensively ukqcd ; mpglue  and the consensus is that the lightest glueball has scalar quantum numbers with the tensor ($`J^{PC}=2^{++}`$) and pseudoscalar glueballs ($`J^{PC}=0^+`$) next in mass. The quenched results have been explored to very small lattice spacings and there is convincing evidence that the continuum limit values have been extracted - see fig. 1. Since the quenched approximation does not reproduce experiment, different ways to set the scale will differ - by $`\pm 10\%`$. Using a conventional scale assignment ($`r_0=0.5`$ fm) gives masses of around 1.6 GeV for the scalar glueball and 2.2 GeV for the tensor.

One signal of great interest would be a glueball with $`J^{PC}`$ not allowed for $`q\overline{q}`$ \- a spin-exotic glueball or oddball \- since it would not mix with $`q\overline{q}`$ states. These states are found MTgl ; ukqcd ; mpglue  to be high lying: considerably above $`2m(0^{++})`$. Thus they are likely to be in a region very difficult to access unambiguously by experiment.

Within the quenched approximation, the glueball states are unmixed with $`q\overline{q},q\overline{q}q\overline{q}`$, etc. Furthermore, the $`q\overline{q}`$ states have degenerate flavour singlet and non-singlet states in the quenched approximation. This gives rise to anomalies in the quenched approximation: for example the scalar meson propagation can have the wrong sign fnal  because the $`\eta \pi `$ intermediate state is mistreated. Once quark loops are allowed in the vacuum then these anomalies are removed. Indeed, for the favour-singlet states of any given $`J^{PC}`$, there will be mixing between the $`s\overline{s}`$ state, the $`u\overline{u}+d\overline{d}`$ state and the glueball. Since the scalar glueball is lightest, we expect the largest effects here and one can explore this by measuring directly the scalar mass eigenstates in a study with $`N_f=2`$ flavours of sea-quark.

Most studies have shown no significant change of the glueball spectrum as dynamical quarks are included bali . However the larger lattice spacing result cmcm176  shows a significant reduction in the lightest scalar mass, as shown in fig. 1. Before concluding that this implies a lower scalar mass in the continuum limit, one needs to check whether an enhanced order $`a^2`$ correction might be present. Studies using the same approach at a finer lattice spacing ht ; cmcm202  do suggest that this large order $`a^2`$ effect is significant, but studies even nearer to the continuum or with improved actions are needed to resolve this fully.

Let us now discuss the mixing of the scalar glueball and scalar mesons. In quenched QCD there is the problem described above concerning scalar mesons. For heavy enough quarks this is unimportant and one can measure the mixing strength on a quenched lattice even though no mixing actually occurs. On a rather coarse lattice ($`a^11.2`$ GeV), two groups have attempted this weinssg ; cmcm176 . Their results expressed as the mixing for two degenerate quarks of mass around the strange quark mass are similar, namely $`E0.36`$ GeV weinssg  and 0.44 GeV cmcm176 . Opinions differ weinssg ; cmcm176  as to whether this large mixing which would shift the glueball mass down by 20% will persist to the continuum limit.

From dynamical fermion studies with $`N_f=2`$, one can determine the flavour singlet and non-singlet mass spectrum. No glueball, as such, can be defined. What we find cmcm176 ; cmcm202  is that the lightest flavour-singlet scalar meson is lighter than the lightest flavour non-singlet. This is in qualitative agreement with the mixing scenario described above.

As well as this mixing of the glueball with $`q\overline{q}`$ states, there will be mixing with $`q\overline{q}q\overline{q}`$ states which will be responsible for the hadronic decays. A first attempt to study this gdecay  in quenched QCD yields an estimated width for decay to two pseudoscalar mesons from the scalar glueball of order 100 MeV. A more realistic study would involve taking account of mixing using gluonic, $`q\overline{q}`$ and $`q\overline{q}q\overline{q}`$ operators in full QCD.

## 4 Hybrid mesons

A hybrid meson has the gluonic degrees of freedom which are excited non-trivially. The most significant consequence of this, experimentally, will be mesons with $`J^{PC}`$ values not allowed in the quark model for a $`q\overline{q}`$ system (such as $`1^+`$).

On the lattice one can easily deal with relatively light quarks (down to the strange quark mass) or extremely heavy quarks (treated as static or using NRQCD). I first discuss hybrid mesons with static heavy quarks where the description can be thought of as an excited colour string. The lattice results liv ; pm ; jkm  can then be presented as potential energy versus quark-antiquark separation $`R`$. The ground state will correspond to the usual interquark potential whereas excited states will have non-trivial representations of the symmetries, e.g. the lightest excited state ($`\mathrm{\Pi }_u`$) has colour flux which is a difference of paths from quark to antiquark of the form $``$.

From the potential corresponding to these excited gluonic states, one can determine the spectrum of hybrid quarkonia using the Schrödinger equation in the Born-Oppenheimer approximation. The $`\mathrm{\Pi }_u`$ symmetry state will produce a degenerate set of eight hybrid mesons of which those with $`J^{PC}=1^+,0^+`$ and $`2^+`$ are spin-exotic and hence will not mix with $`Q\overline{Q}`$ states. They thus form a very attractive goal for experimental searches for hybrid mesons.

Within the quenched approximation, the lattice evidence for $`b\overline{b}`$ quarks points to a lightest hybrid spin exotic meson H with $`J^{PC}=1^+`$ at an energy given by hf8  $`(m_Hm_{2S})r_0=1.9\pm 0.1`$. This has been checked with $`N_f=2`$ studies cppacs ; hdecay  which find similar mass ratios $`(m_Hm_{1S})/(m_{1P}m_{1S})`$ but a larger splitting in terms of $`r_0`$, namely cppacs  $`(m_Hm_{2S})r_0=2.4\pm 0.2`$. Using the experimental mass of the $`\mathrm{{\rm Y}}(2S)`$, these results imply that the lightest spin exotic hybrid is in the mass range 10.7 to 10.9 GeV. Above this energy there will be many more hybrid states, many of which will be spin exotic.

Within this static quark framework, one can explore the decay mechanisms. One special feature is that the symmetries of the quark and colour fields about the static quarks must be preserved exactly in decay. This has the consequence that the decay from a $`\mathrm{\Pi }_u`$ hybrid state to the open-$`b`$ mesons ($`B\overline{B},B^{}\overline{B},B\overline{B^{}},B^{}\overline{B^{}}`$) will be forbidden hdecay  if the light quarks in the $`B`$ and $`B^{}`$ mesons are in an S-wave relative to the heavy quark. The decay to $`B^{}`$-mesons with light quarks in a P-wave is allowed by symmetry but not energetically. The only allowed decays are to $`\chi _b+M`$ where $`M`$ is a light quark-antiquark meson in a flavour singlet. Lattice estimates hdecay  of these transitions have been made and the dominant mode (with a width of around 100 MeV) is found to be with $`M`$ as a scalar meson. This decay analysis does not take into account heavy quark motion or spin-flip and these effects will be significantly more important for charm quarks than for $`b`$-quarks.

Several lattice groups livhyb ; milc ; sesamhyb ; ml  have studied hybrid spectra for light quarks and find the lightest spin-exotic hybrid to have $`J^{PC}=1^+`$ and mass (for $`s`$-quarks) of around 2GeV. The corresponding light-quark ($`u`$, $`d`$) state would be around 120 MeV lighter. The light quark results have also been milc  extrapolated to charm quarks and masses near 4.4GeV are found for the corresponding state. A recent study has confirmed this value but with a very long extrapolation ml . These mass estimates can be compared to naive estimates pm  of the spin-exotic charm state mass of 4.0 GeV from the static quark approach which will have an uncontrolled systematic error.

It is not easy to reconcile these lattice results with experimental indications expt  for resonances at 1.4 GeV and 1.6 GeV, especially the lower mass value. Mixing with $`q\overline{q}q\overline{q}`$ states such as $`\eta \pi `$ is not included for realistic quark masses in the lattice calculations. This can be interpreted, dependent on one’s viewpoint, as either that the lattice calculations are incomplete or as an indication that the experimental states may have an important meson-meson component in them.

An attractive prospect to study hybrid mesons is from $`p\overline{p}`$ annihilation. Spin-exotic hybrids, which provide the best signal, will only be produced in association with other hadrons. The decay channel to $`\chi _c+\pi +\pi `$ with the two $`\pi `$ mesons in an S-wave is a promising detection channel. Lattice studies give first estimates of the masses and decay widths but more work needs to be done to constrain the systematic errors on these estimates.


text:
# A High-Resolution Study of 𝜼 Carinae’s Outer Ejecta

## 1. Morphology, kinematics and the soft X-ray emission

The outer part of the nebula around $`\eta `$ Carinae contains a countless number of knots, bullets and filaments and manifests the so-called outer ejecta (diameter 60′′ or 0.67 pc). A kinematic analysis using high-resolution echelle spectra showed that the outer ejecta expand with velocities between $`600`$ km s-1 and +600  km s-1 on average (e.g. Meaburn et al. 1996, Weis 2001a,b). Various filaments nevertheless show much higher radial velocities—e.g. $``$ 2000  km s-1(Weis 2001a,b). It is still unclear what triggered the outburst in 1843 and which mechanism not only formed $`\eta `$ Carinae’s nebula but also the amazingly high expansion velocities. With velocities that high, X-ray emission from $`\eta `$ Carinae is expected and was indeed detected (e.g. Cheblowski et al. 1984, Corcoran et al. 1997). With ROSAT and CHANDRA we are now able to separate between a harder, nearly point like emission from the central source and an extended softer emission from $`\eta `$ Carinae’s outer ejecta which is roughly hook shaped (e.g. Weis et al. 2001). Interestingly, an overlay of the optical image and the X-ray emission shows only a few correlations between the optically emitting gas and the hot X-ray gas, see Fig. 1. A much stronger conformance was found in comparison with the radial expansion velocities. The expansion velocities are derived from our optical echelle spectra (FWHM $``$ 14  km s-1) and overplotted in Fig. 1 with negative velocities underlined. Areas with higher X-ray emission show in general higher expansion velocities. The faster the gas is expanding the more intense is the X-ray emission. The extended soft X-ray emission from $`\eta `$ Carinae and especially the morphology of the X-ray nebula can therefore be explained by faster moving filaments which form X-ray emitting shocks. The temperature of the gas (using thermal equilibrium models) is $``$ 0.65 keV (Weis et al. 2002, in prep.) indicating post-shock velocities of around 750  km s-1 in very good agreement with the detected expansion velocities of the bulk of the gas.

## References

Chlebowski, T., Seward, F. D., Swank, J., Szymkowiak, A., 1984, ApJ, 281, 665

Corcoran, M. F., Ishibashi, K., Swank, J. H., et al. 1997, Nature, 390, 587

Meaburn, J., Boumis, P., Walsh, J. R., et al. 1996, MNRAS, 282, 1313

Weis, K., 2001a, in $`\eta `$ Carinae and other mysterious stars, ed. T. Gull, S. Johanson, & K. Davidson, ASP Conf. 242, 129

Weis, K., 2001b, in Reviews in Modern Astronomy 14, Springer-Verlag, ed.: Schielicke, R.E, ISBN 3-9805176-4-0, 261

Weis, K., Duschl, W.J., Bomans, D.J., 2001, A&A, 367, 566


text:
# Monochromatic path crossing exponents and graph connectivity in 2D percolation

## Acknowledgements

We thank S. Kirkpatrick and J. Vannimenus for useful discussions.


text:
# Entropy production rates of bistochastic strictly contractive quantum channels on a matrix algebra

## References


text:
# Dynamical nature of the nuclear pseudospin and its isospin asymmetry

## 1 Introduction

Pseudospin symmetry was introduced in the late 60’s to account for a degeneracy in single-particle nucleon levels of certain nuclei, with quantum numbers ($`n`$, $`\mathrm{}`$, $`j=\mathrm{}+1/2`$) and ($`n1`$, $`\mathrm{}+2`$, $`j=\mathrm{}+3/2`$) . These levels have the same quantum number $`\stackrel{~}{\mathrm{}}=\mathrm{}+1`$, called “pseudo” orbital angular momentum, and a “pseudo” spin quantum number, $`\stackrel{~}{s}=1/2`$. In general, for a $`n\mathrm{}j`$ state, we have $`\stackrel{~}{\mathrm{}}=\mathrm{}\pm 1`$ if $`j=\mathrm{}\pm 1/2`$. Pseudospin symmetry is exact when doublets with $`j=\stackrel{~}{\mathrm{}}\pm \stackrel{~}{s}`$ are degenerate. Since then the origin of pseudospin symmetry in nuclei has been a subject of intense theoretical research.

Recently, Ginocchio  uncovered the relativistic character of the symmetry, noting that the pseudo-orbital angular momentum is the orbital angular momentum of the lower component of the Dirac spinor. This is an exact SU(2) symmetry for the Dirac Hamiltonian  with an attractive scalar potential, $`S`$, and a repulsive vector potential $`V`$, when they have the same magnitude, i.e., $`\mathrm{\Sigma }=S+V=0`$. This would explain why there are quasi-degenerate pseudospin doublets in nuclei, since in relativistic mean field (RMF) models nuclear saturation arises from an extensive cancellation between a large attractive scalar potential and a large repulsive vector potential. However, since $`\mathrm{\Sigma }`$ acts as a binding potential, one cannot set $`\mathrm{\Sigma }=0`$ in the nucleus, because then there would not exist any bound states. Other authors have studied the role of pseudospin-orbit coupling in pseudospin symmetry .

## 2 Dynamical pseudospin symmetry and isospin asymmetry

Pseudospin symmetry is related to an invariance under a SU(2) transformation of the Dirac Hamiltonian with central scalar $`S`$ and vector $`V`$ potentials

$`H=𝜶\mathbf{}𝒑+\beta (m+S)+V,`$ (1)
whose generators are  $`S_i=s_i\frac{1}{2}(1\beta )+\stackrel{~}{s}_i\frac{1}{2}(1+\beta )`$ where $`\stackrel{~}{s}_i=\frac{𝝈𝒑}{p}s_i\frac{𝝈𝒑}{p}`$, $`s_i=\sigma _i/2`$. The commutator of $`S_i`$ with (1) is zero if $`\mathrm{\Sigma }=0`$ or $`\frac{d\mathrm{\Sigma }}{dr}=0`$ .

The Dirac equation $`H\mathrm{\Psi }=ϵ\mathrm{\Psi }`$ can be written as a pair of second-order differential equations for the upper and lower components. Setting $`\mathrm{\Delta }=VS`$, $`V`$ and $`S`$ being radial functions, and defining $`ϵ=E+m`$, we have

$`p^2\mathrm{\Psi }_+`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }^{}}{E+2m\mathrm{\Delta }}}\left({\displaystyle \frac{}{r}}{\displaystyle \frac{1}{r}}𝝈𝑳\right)\mathrm{\Psi }_++(E+2m\mathrm{\Delta })(E\mathrm{\Sigma })\mathrm{\Psi }_+`$ (2)
$`p^2\mathrm{\Psi }_{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Sigma }^{}}{E\mathrm{\Sigma }}}\left({\displaystyle \frac{}{r}}{\displaystyle \frac{1}{r}}𝝈𝑳\right)\mathrm{\Psi }_{}+(E+2m\mathrm{\Delta })(E\mathrm{\Sigma })\mathrm{\Psi }_{},`$ (3)
where $`\mathrm{\Psi }_\pm =[(1\pm \beta )/2]\mathrm{\Psi }`$ are the upper and lower components and the primes denote derivatives with respect to $`r`$. The $`𝝈𝑳`$ term in (3) is the pseudospin-orbit term.

We solve the Dirac equation with the Hamiltonian (1) with scalar and vector potentials with Woods-Saxon shape: $`U(r)=U_0/\left(1+\mathrm{exp}[(rR)/a]\right)`$. We fitted the parameters of this potential to the neutron spectra of 208Pb, obtaining the values $`R=7`$ fm, $`\mathrm{\Delta }_0=650`$ MeV, $`\mathrm{\Sigma }_0=66`$ MeV and $`a=0.6`$ fm. We then varied the diffuseness, radius and $`\mathrm{\Sigma }_0`$ separetely, keeping all the other parameters fixed and showed that the splittings of the pseudospin doublets vary notably, namely they decrease as $`|\mathrm{\Sigma }_0|`$ decreases and diffuseness increases, sometimes reversing sign. The details can be found in ref. .

The effect of the pseudospin-orbit term on pseudospin energy splittings can be assessed by writing (3) as a Schroedinger-like equation, dividing it by an energy- and $`r`$-dependent effective mass $`m^{}=(E+2m\mathrm{\Delta })/2`$ and taking the expectation value of each term relative to $`\mathrm{\Psi }_{}`$:

$$\frac{p^2}{2m^{}}+V_{\mathrm{PSO}}+V_\mathrm{D}+\mathrm{\Sigma }=E,$$
(4)
where

$`{\displaystyle \frac{p^2}{2m^{}}}`$ $`=`$ $`{\displaystyle \mathrm{\Psi }_{}^{}\frac{p^2}{2m^{}}\mathrm{\Psi }_{}\mathrm{d}^3𝒓/\mathrm{\Psi }_{}^{}\mathrm{\Psi }_{}\mathrm{d}^3𝒓},`$ (5)
$`V_{\mathrm{PSO}}`$ $`=`$ $`\mathrm{P}{\displaystyle \mathrm{\Psi }_{}^{}\frac{1}{2m^{}}\frac{\mathrm{\Sigma }^{}}{E\mathrm{\Sigma }}\frac{1}{r}𝝈𝑳\mathrm{\Psi }_{}\mathrm{d}^3𝒓/\mathrm{\Psi }_{}^{}\mathrm{\Psi }_{}\mathrm{d}^3𝒓}`$ (6)
$`V_\mathrm{D}`$ $`=`$ $`\mathrm{P}{\displaystyle \mathrm{\Psi }_{}^{}\frac{1}{2m^{}}\frac{\mathrm{\Sigma }^{}}{E\mathrm{\Sigma }}\frac{\mathrm{\Psi }_{}}{r}\mathrm{d}^3𝒓/\mathrm{\Psi }_{}^{}\mathrm{\Psi }_{}\mathrm{d}^3𝒓},`$ (7)
$`\mathrm{\Sigma }`$ $`=`$ $`{\displaystyle \mathrm{\Psi }_{}^{}\mathrm{\Sigma }\mathrm{\Psi }_{}\mathrm{d}^3𝒓/\mathrm{\Psi }_{}^{}\mathrm{\Psi }_{}\mathrm{d}^3𝒓},`$ (8)
where ‘P’ denotes the principal value of the integral. These terms can be identified, respectively, as a kinetic term, a pseudospin-orbit term, a potential term related to what is sometimes called Darwin term and the mean value of the $`\mathrm{\Sigma }`$ potential with respect to the lower component, $`\mathrm{\Psi }_{}`$.

In Fig. 1 are plotted the differences of the terms in Eq. (4) between each member of the three topmost pseudospin partners in 208Pb, $`(1i_{11/2},2g_{9/2})`$, $`(2f_{5/2},3p_{3/2})`$ and $`(1h_{9/2},2f_{7/2})`$. One sees that the contribution of $`V_{\mathrm{PSO}}`$ for the pseudospin energy splittings is larger than the splittings themselves and has the opposite sign of the kinetic and $`\mathrm{\Sigma }`$ terms. For all doublets there is a significant cancellation of $`V_{\mathrm{PSO}}`$ with the kinetic and $`\mathrm{\Sigma }`$ terms. A similar analysis was made in ref. , concluding that the pseudospin-orbit potential is non-perturbative. We also found that there is a clear correlation between $`V_{\mathrm{PSO}}`$ and the pseudospin energy splitting when the diffusivity and the depth $`\mathrm{\Sigma }_0`$ are varied .

All these findings point to the conclusion that pseudospin symmetry has a dynamical nature, since it is very much dependent on the shape of the $`\mathrm{\Sigma }`$ nuclear mean field potential and results from a strong cancellation among the several contributing terms to the energy splittings.

This pseudospin dependence on the shape of the $`\mathrm{\Sigma }`$ potential allows us to explain why the proton and neutron levels have different pseudospin splittings .

In RMF models, the $`\rho `$ meson interaction modifies the vector potential to $`V=V_\omega +V_\rho =V_\omega \pm \frac{g_\rho }{2}\rho _0`$, where the plus sign refers to protons and the minus sign refers to neutrons, and $`\rho _0`$ is the time component of the $`\rho `$ field which is proportional to $`ZN`$. The magnitude of the $`\rho `$ potential is around 4-8 MeV. For heavy nuclei, $`\rho _0<0`$ and therefore $`V_p<V_n`$ which implies $`|\mathrm{\Sigma }_n|<|\mathrm{\Sigma }_p|`$. This means that, for neutron-rich nuclei, the pseudospin symmetry is favored for neutrons. Notice, however, that the Coulomb potential, which has the opposite sign of $`V_\rho `$ for protons, also affects proton spectra. Both of these effects can be seen in Figure 2, where is plotted the effect of $`V_\rho `$ and the Coulomb potential on the energy splitting of the $`(2d_{5/2},1g_{7/2})`$ pseudospin partners along the Sn isotopic chain from $`A=130`$ to $`A=170`$. We see that the potential $`V_\rho `$ gives the main contribution for the observed isospin asymetry.

Acknowledgments

We acknowledge financial support from FCT (POCTI), Portugal, and from CNPq/ICCTI Brazilian-Portuguese scientific exchange program.


text:
# A Warp in the LMC Disk?

## 1 Introduction

Among the many virtues of the Large Magellanic Cloud (LMC), one of the most important is that in spite of its proximity, for most purposes its stars may be considered to lie at a single distance. As a result, the LMC is frequently used in studies of stellar evolution, is of primary importance in the calibration of the distance scale, and is the target of several projects aimed at detecting microlensing events in the Galactic halo (e.g. Alcock et al. 2000, Lasserre et al. 2000, Udalski et al. 1999, Drake et al. 2002). The neighboring Small Magellanic Cloud (SMC) has never reached the same level of popularity, for the interesting reason that it is stretched by several kpc along the line of sight (e.g. Caldwell & Coulson 1986, Hatzidimitriou & Hawkins 1989), the probable result of tidal interactions with the LMC and the Milky Way.

Probing deeper, we find that the stars of the LMC are not, after all, quite at the same distance, because the LMC disk is tilted with respect to the line of sight. The LMC’s tilt has most recently been measured by van der Marel & Cioni (2001), who found its value to be $``$35° using asymptotic giant branch and red giant branch tip (TRGB) stars observed in the 2MASS survey (Skrutskie 1998). The effect of the tilt is to produce a modulation of the apparent luminosity of these tracers with position angle in the LMC’s disk, with an amplitude of $``$0.2 magnitude peak-to-trough. Moreover, upon deprojecting the 2MASS data, van der Marel (2001) found that the LMC disk is elliptical and has a nonuniform surface density distribution, indicative of tidal forces acting on the LMC. Clearly, the LMC is not slipping through its interactions with the Milky Way and the SMC completely unscathed.

We are thus led to ask, what further evidence for the LMC-SMC-Milky Way interaction may we find hidden in the LMC’s geometry? The question is important for chronicling the interaction history through numerical models (e.g. Gardiner & Noguchi 1996), and may have implications for interpreting the microlensing results (e.g. Zaritsky & Lin 1997, Zhao & Evans 2000). In this paper, we report the results of our use of the apparent luminosity of “red clump” stars to measure the shape of the LMC disk. The red clump, which is composed of the younger, more metal-rich analogues to the core helium-burning horizontal branch stars found in globular clusters, has been widely discussed in the literature for use as a distance indicator (e.g. Girardi & Salaris 2001, GS01; Cole 1998). While the red clump remains controversial as an absolute distance indicator, by using it to measure relative distances to LMC fields we avoid the disagreement over its zero point. Red clump stars are more numerous than TRGB stars by a factor of $``$100 and have tightly defined colors and magnitudes, making them easily identifiable in color-magnitude diagrams (CMDs). In the LMC, they have a typical surface density of 3.5$`\times 10^4`$ stars per square degree, allowing us to search for structure with high spatial frequency. Because they have well-defined mean colors in addition to magnitudes, we can measure reddenings with the same tracers used to measure distances, avoiding the population-dependent reddening effects discussed by Zaritsky (1999). The errors associated with using the red clump as a distance indicator are mainly systematic ones. We discuss their possible effect on our results in section 4.

## 2 Observations

We observed 50 fields in the LMC with the CTIO 0.9-m telescope and SITe 2K$`\times `$2K CCD#3 camera during the nights 2001 November 18–24 (Fig. 1). We used CTIO’s Johnson $`V`$ and Kron-Cousins $`I`$ filters. The fields were randomly selected from a 6°$`\times `$6°area of the LMC, with the constraint that they should lie far from regions of recent star formation, where we expect extinction by interstellar dust to be high. We used gain setting 2, which provides gain of $``$1.5 $`e^{}`$/ADU and read noise of $``$3.5 $`e^{}`$. All of the nights were photometric, while the seeing was typically 1″–1$`\stackrel{}{\mathrm{.}}`$5. Starting with our first field, our procedure was to take one pair of $`V`$ and $`I`$ exposures in one field before moving on to the next. Once all fields had been observed once in $`V`$ and $`I`$, we repeated the sequence twice. With 250 sec of exposure in $`V`$ and 200 sec in $`I`$, we were able to achieve S/N$`100`$ at $`V=19,I=18`$ while observing $`2025`$ fields per night. In all, we obtained 3–4 pairs of $`VI`$ images per field, each observed on different nights.

We observed Landolt (1992) photometric standards at several times during each night for calculation of zero points, atmospheric extinctions, and color corrections. We also observed two standard star fields located near the position of the LMC, for which Arlo Landolt kindly provided standard magnitudes. These fields were used to check for different atmospheric extinction at the LMC’s position compared to the celestial equator, for which we found no evidence.

A byproduct of the fast quad amplifier readout used with the SITe 2K$`\times `$2K CCD is crosstalk between each amplifier and the other three, which produces faint ghosts in the images. Our first image reduction step was to remove these ghosts by subtracting from each CCD section the proper fraction of the neighboring sections. Following this crosstalk correction, we performed the standard image reduction steps of overscan subtraction, trimming, bias subtraction, and flat fielding with twilight sky flats, all done with the QUADPROC package in IRAF.

## 3 Analysis

### 3.1 Photometry

We performed photometry on each of the 304 individual reduced frames with DAOPHOT (Stetson 1987) and the accompanying profile-fitting program ALLSTAR. Because finding isolated stars in the crowded LMC fields was difficult, we selected a large number ($``$200) of stars per frame to define the point-spread function (PSF), allowing an average, quadratically variable PSF to be constructed. We corrected the PSF magnitudes to an aperture of diameter 15″ by performing aperture photometry on the PSF stars after subtracting all other stars, while for some stars allowing the magnitude to be extrapolated based on smaller apertures using the aperture growth curve measured using DAOGROW (Stetson 1990).

### 3.2 Photometric calibration

We measured only aperture magnitudes for the standard star frames, again using DAOGROW to extrapolate the measurements of some of the stars to a 15″ aperture. We then calculated nightly zero points, atmospheric extinction coefficients, and color terms. After averaging the seven color term measurements, we recalculated the zero points and atmospheric extinctions while holding the color terms fixed. We applied these calibration equations to our LMC frames using the PSF stars as local standards. Table 1 shows the calibration coefficients of the transformation equations:

$$v=V+A_0+A_1(VI)+A_2X$$

$$i=I+B_0+B_1(VI)+B_2X$$

where uppercase refers to standard magnitudes, lowercase to instrumental magnitudes, and $`X`$ is the airmass.

### 3.3 Measurement of red clump stars

Each frame contains several thousand red clump stars, easily identified by their location in the CMD. After selecting stars within a box encompassing the red clump (Fig. 2), we fit a Gaussian profile combined with a second order polynomial to the $`VI`$ and $`I`$ histograms within the box, taking the Gaussian means as the mean color and magnitude of the red clump stars. We used bin sizes corresponding approximately to the precision which we can measure the means, given the intrinsic $`VI`$ and $`I`$ dispersions of the histograms and the number of red clump stars per frame. We filtered the histograms with a 20–40 bin-wide Savitzky-Golay filter (Press et al. 1992) before fitting, as we found that our fitting routine behaved more reliably with this step included. However, the $`\chi ^2`$ values and likelihoods of the fits were calculated using the unfiltered histograms. On average, we found the fits to the $`VI`$ histograms produced $`\chi ^2/\nu `$=1.3 with $``$20% chance of measuring higher $`\chi ^2`$ purely by chance; for the $`I`$ histograms the numbers were $`\chi ^2/\nu `$=1.04 with $``$35% probability for higher $`\chi ^2`$.

Next, we averaged the 3–4 measurements per field to produce an average red clump $`VI`$ color and $`I`$ magnitude for each field. Assuming an intrinsic red clump color of $`VI`$=0.92, we calculated $`E(VI)`$ interstellar reddening values for each field and converted these to $`I`$-band interstellar extinctions through the equation $`A_I=1.4E(VI)`$ (Schlegel, Finkbeiner, & Davis 1998). We selected the intrinsic color so as to produce a median reddening equal to that measured by Schlegel et al. (1998) for the LMC, but note that our choice has no effect on our results since we are using the red clump only to measure relative distances.

Our formal errors in measuring the mean color and magnitude of the red clumps are small, $`<<`$1%. We thus expect our distance errors to be dominated by the photometric calibration (1%) and by error in assigning the correct reddening through the red clump color (2%). We expect that variations in the age and metallicity of the red clump from field to field will be $``$2%; however, we do not include this source in our budget of random distance errors, as it is likely to be correlated with position in the LMC. We discuss the sources of error in depth in section 4.

### 3.4 The LMC’s geometry

For the purposes of this section, we assume that the variation in the reddening-corrected, mean $`I`$ magnitudes of the red clump stars in our fields is produced purely by distance effects, so that we can study the LMC’s geometry. First, we converted the celestial coordinates of our fields into planar coordinates following Wesselink (1959). Second, we transformed the observed $`\mathrm{\Delta }I=I^{RC}I_{}^{RC}`$ into distances relative to the center of the LMC, which we take to be 05h19m38$`\stackrel{\mathrm{s}}{\mathrm{.}}`$0 $``$69°27′5$`\stackrel{}{\mathrm{.}}`$2 (2000.0, precessed from 1950.0; de Vaucouleurs & Freeman 1973). We note that our assumed origin has no effect on our study of the LMC’s shape. Finally, we corrected for the small geometric projection effect produced by our fixed vantage point.

A weighted least-squares fit of a plane to all of the data yields a disk inclination angle $`i=24\mathrm{°}\pm 2\mathrm{°}`$ and position angle of the line of nodes $`\theta =149\mathrm{°}\pm 6\mathrm{°}`$, where we have used the convention of $`i=0\mathrm{°}`$ for a face-on disk and measured the position angle as the number of degrees east with respect to north. This measurement is in good agreement with e.g. Caldwell & Coulson (1986), who found $`i=28\stackrel{}{\mathrm{.}}6\pm 5\stackrel{}{\mathrm{.}}9`$ and $`\theta =142\stackrel{}{\mathrm{.}}4\pm 7\stackrel{}{\mathrm{.}}7`$ from their analysis of individual Cepheid distances. However, there is a possible problem with the fit, as $`\chi ^2`$ per degree of freedom is equal to 1.4, with only a 2% chance of attaining $`\chi ^2`$ as high as ours through chance alone.

Figure 3 (right panel) shows the LMC disk as viewed edge-on along our measured line of nodes. The tilt of the disk is strikingly well defined, particularly when compared to Figure 6 of Caldwell & Coulson (1986). However, the fields at the southwestern edge (open circles) fall systematically low compared to the fitted plane, by $``$0.1 magnitudes in $`I`$ or 2.5 kpc in distance. If we remove these 15 fields (labelled by diamonds in the left panel of Figure 3) from the fit, we find instead $`i=35\stackrel{}{\mathrm{.}}8\pm 2\stackrel{}{\mathrm{.}}4`$ and $`\theta =145\mathrm{°}\pm 4\mathrm{°}`$. This inclination value is in excellent agreement with the value $`i=34\stackrel{}{\mathrm{.}}7\pm 6\stackrel{}{\mathrm{.}}2`$ measured by van der Marel & Cioni (2001), although we disagree at the $`2\sigma `$ level with their position angle $`\theta =122\stackrel{}{\mathrm{.}}5\pm 8\stackrel{}{\mathrm{.}}3`$. The statistical likelihood of the fit is now excellent, with a 96% chance of producing $`\chi ^2`$ values higher than our measured 0.62 per degree of freedom. Indeed, the high likelihood of the fit suggests that we have overestimated the size of our error bars.

Fitting a plane to only the 15 southwestern fields which we removed above, we find $`i=14\mathrm{°}\pm 9\mathrm{°}`$ and $`\theta =280\mathrm{°}\pm 47\mathrm{°}`$, with $`\chi ^2/\nu =0.5`$ and 94% probability of producing higher $`\chi ^2`$ by chance. Thus, our data imply that the LMC disk is warped and that its line of nodes is twisted at the southwestern edge.

## 4 Discussion

As shown by fitting planes to both the 15 southwestern and 35 remaining fields, our analysis suggests that the disk of the LMC is warped and twisted. This result can perhaps explain the different tilts measured by Caldwell & Coulson (1986) and van der Marel & Cioni (2001) as being caused by differing contributions of the warp to their datasets. The area sampled by our observations is similar to that studied by Caldwell & Coulson (1986), so that we expect similar results if we include all of our fields in the fit. Moreover, van der Marel & Cioni (2001) find that the inclination of the LMC drops by $``$10 degrees outside a radius of $``$4$`\stackrel{}{\mathrm{.}}`$5 from the center. Although our study suggests that the LMC’s warp begins at a radius $``$3°from the center, van der Marel & Cioni (2001) surveyed a much larger area than we have, making a direct comparison with their study difficult.

Our finding that the LMC is warped depends, however, on our assumption that distance is the dominant effect producing the observed variation in red clump luminosity. In the following sections, we discuss the possible sources of systematic error.

### 4.1 Uncertainty in the calibration

The stability of the derived photometric zero points and atmospheric extinction terms suggests that the photometric calibration is accurate to $``$1% (Table 1). A further check against systematic errors is provided by the fact that we observed every field on three different nights. In Figure 4, we plot $`VI`$ and $`I`$ of the red clump for the 152 pairs of observations, having subtracted the mean clump $`VI`$ and $`I`$ for each field. The nightly drift and scatter is smaller than 1%, ruling out calibration uncertainty as the cause of the apparent warp.

### 4.2 Reddening

If we have overestimated the reddening $`E(VI)`$ in the 15 southwestern fields by $``$0.07 magnitudes compared to the other fields, then this could explain the effect we see in Figure 3. Figure 5, where we plot $`A_I`$ as derived from the color of the red clump versus the departure of the fields in $`I`$ from the plane with $`i=35\mathrm{°}`$, is cause for suspicion. The southwestern fields have red clumps with systematically higher interstellar extinction compared to the other fields. Moreover, the magnitude deviation from the plane is correlated with $`A_I`$, such that if we have overestimated the reddening in the southwestern fields by more than a factor of two, then reddening could explain the supposed warp.

We checked whether we have miscalculated the reddening by looking for differences in other evolutionary features in the color-magnitude diagrams of the fields. If the reddening values are wrong, the error should manifest itself as a difference in color between the red clump and the red giant branch and main sequence in the southwestern fields. After correcting for the reddening and red clump luminosity differences, we merged the photometry for the southwestern fields and the remaining fields separately. Figure 6 shows the difference of the two-dimensional histograms of the CMDs of the merged southwestern and other fields, after applying various shifts in $`VI`$ to the southwestern fields. The overlay of the colors of the main sequence and red giant branch allow a $``$0.01–0.02 magnitude error in $`E(VI)`$, but not a difference as large as 0.07 magnitudes. We therefore rule out large reddening errors as the cause of the warp.

If we have properly accounted for interstellar correction, then why is there a correlation between the magnitude deviation from the LMC plane and $`A_I`$? Figure 7 offers a plausible explanation. The COBE/DIRBE–IRAS/ISSA dust map (Schlegel et al. 1998) shows that the southwestern corner of the LMC is associated with a region of higher diffuse extinction, presumably due to Galactic foreground dust, having $`A_I0.3`$. For comparison, we measured $`A_I0.25`$ at the LMC’s southwest edge.

### 4.3 Age and metallicity

GS01 and Cole (1998) have discussed the strong effect that age and metallicity have on the luminosity of red clump stars. In short, metal-poor clumps are brighter than metal-rich ones, and younger clumps brighter than older ones. Assuming a mean red clump age of $``$4 Gyr, appropriate for the LMC (GS01), and a mean metallicity of \[Fe/H\]$`0.6`$ (Pagel & Tautvaisiene 1998), $`M_I`$ of the red clump may be made brighter by 0.1 magnitude by either lowering \[Fe/H\] by $``$1 dex or by lowering the age by $``$2 Gyr; $`M_I`$ could be made fainter by 0.1 magnitude by raising \[Fe/H\] by 0.3 dex or by increasing the age by $``$3 Gyr (GS01). Thus, we need to consider that the supposed warp may be produced by the southwestern fields being more metal-poor and/or younger than the rest of the LMC disk.

There are a number of arguments against this hypothesis. First, the appearance of the CMDs appear to rule out large bulk differences in age or metallicity between the southwestern and the remaining fields. Figure 8 shows the combined CMDs plotted against Girardi et al. (2000) isochrones with ages of 2 and 4 Gyr and \[Fe/H\] equal to $`0.6`$ and $`1.6`$. The similar slopes of the red giant branches in the southwestern fields compared to the other fields suggest that the ages and metallicities of the red clump stars are also similar.

Second, the variation due to age and metallicity of the $`I`$ luminosity of the red clump is predicted to be small across the LMC. GS01, using the LMC Bar and outer disk star formation histories derived by Holtzman et al. (1999) coupled with the age-metallicity relation predicted by Pagel & Tautvaisiene (1998), calculated that the Bar red clump should be only 0.01 $`I`$ magnitudes brighter than the outer disk red clump. We derive a similar result by feeding the star formation histories derived by Smecker-Hane et al. (2002) for the LMC Bar and disk through the GS01 models. We find that the LMC Bar should have a red clump brighter than the disk by $``$0.03 magnitudes, despite the Bar having a $``$2 Gyr younger age. For comparison, Smecker-Hane et al. (2002) find that the red clump in their Bar field is 0.02 magnitudes brighter in $`I`$ than in their disk field. After correcting for the 0.04 magnitude larger $`A_I`$ in their Bar field and for our solution of the LMC’s tilt, which places their Bar field 0.02 magnitudes closer than their disk field, we find that their Bar red clump is 0.04 magnitudes more luminous than their disk red clump. This measurement is in good agreement with the theoretical prediction. Thus, both theory and observations show that the mean properties of field red clump stars are less sensitive to age and metallicity effects than are single-age, single-metallicity populations. This is to be expected from the composite nature of stars in the field.

Third, differential rotation in the LMC shreds unbound coherent structure on the orbital timescale (270 Myr at $`R`$=3 kpc, adopting $`v_{\mathrm{rot}}=70`$ km s-1; Alves & Nelson 2000), while the bulk of the red clump stars in the LMC are at least $``$4 Gyr old (GS01). We thus do not expect that an area of the LMC $``$2 kpc in size would retain a record of its locally produced age and metallicity distribution, distinct from fields at similar radii. A warp, on the other hand, would be a non-equilibrium, recently formed structure not yet mixed into the surroundings.

Finally, explaining the supposed warp through population effects would require an age or metallicity gradient that is inconsistent with observed trends. The general trend in the LMC is for the younger clusters to lie preferentially towards the inner regions (Bica et al. 1996), in the oppposite sense from that required here. The required metallicity gradient ($`1.0`$ dex over a stretch of $``$2 kpc), while operating in the right sense, is much steeper than the weak radial gradient observed (Kontizas et al. 1993).

In summary, we think that the simplest explanation for the unusually bright red clump luminosities seen in our southwestern fields is that the LMC disk is warped, rather than that the southwest of the LMC is more metal-poor or younger than the rest of the disk. However, a detailed solution of the star formation and chemical enrichment histories implied by our fields is needed to answer the question definitively.

## 5 Conclusions and future work

We have studied the shape of the LMC by measuring the variation in luminosity of red clump stars across 50 13′$`\times `$13′ fields distributed over a 6°$`\times `$6°area. Our observations were taken with the CTIO 0.9-m telescope CCD camera using $`V`$ and $`I`$ filters. We clearly detect the tilt of the LMC disk through a decrease in brightness of the red clump stars along the LMC’s NE-SW axis. Our measured tilt of $`35\stackrel{}{\mathrm{.}}8\pm 2\stackrel{}{\mathrm{.}}4`$ is in excellent agreement with the recent measurement of van der Marel & Cioni (2001), as long as we exclude 15 fields in the southwest which appear to lie out of the plane of the disk. Including all of the fields, we measure a lower inclination of $`24\mathrm{°}\pm 2\mathrm{°}`$. This value is in good agreement with that measured by Caldwell & Coulson (1986); however, the fitted plane is of lower statistical significance.

We believe that we have detected a warp in the LMC disk which causes a region $``$2 kpc wide in the southwest to lie out of the plane by $``$2.5 kpc. If it exists, this out-of-plane feature is likely a non-virialized structure. How and when was it produced? The location of the warp points a finger at the SMC. According to the model by Gardiner & Noguchi (1996), the LMC and SMC endured a close encounter $``$200 Myr ago which drew out the material that new occupies the inter-Cloud region. It is possible that this interaction also altered the shape of the LMC disk. However, the position of the warp also aligns it with the major axis of the LMC’s elliptical disk (van der Marel 2001), the shape of which van der Marel attributes to the tidal influence of the Milky Way. Perhaps our suggested warp is an additional response to the Milky Way’s tidal force? Weinberg (1998, 2000) studied the mutual influence of the LMC and Milky Way through numerical simulations. Although Weinberg found that the long-term effect on the LMC’s structure is simply heating of the disk, transient structures might also be produced.

Further exploration of the structure of the LMC could determine whether the warp is reflected in northeastern fields more distant than those observed here. The determination of the star formation and chemical enrichment history in the LMC’s southwest would establish with certainty the contribution of age and metallicity effects to the observed red clump luminosities. In this study we avoided the LMC’s inner regions, including the Bar, as these are heavily crowded. Measurement of the LMC’s shape and thickness in these regions is of strong interest for understanding the contribution of LMC self-lensing to the observed microlensing events (Zhao & Evans 2000). While our results bear no consequence for the microlensing studies, they do add to the growing body of evidence suggesting that the LMC is not a simple flattened disk. A study of the red clump in the LMC Bar, combined with the continued monitoring of the LMC for microlensing events by the SuperMACHO project (Drake et al. 2002), could provide very interesting results.

This work was carried out as part of the 2002 Research Experiences for Undergraduates (REU) program at CTIO, which is supported by the National Science Foundation (NSF). CTIO is part of the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under cooperative agreement with the NSF. We acknowledge useful comments from the anonymous referee.