text:
# Threshold effects in excited charmed baryon decays

## Abstract

Motivated by recent results on charmed baryons from CLEO and FOCUS, we reexamine the couplings of the orbitally excited charmed baryons. Due to its proximity to the $`\mathrm{\Sigma }_c\pi `$ threshold, the strong decays of the $`\mathrm{\Lambda }_c^+(2593)`$ are sensitive to finite width effects. This distorts the shape of the invariant mass spectrum in $`\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi ^+\pi ^{}`$ from a simple Breit-Wigner resonance, which has implications for the experimental extraction of the $`\mathrm{\Lambda }_c^+(2593)`$ mass and couplings. We perform a fit to unpublished CLEO data which gives $`M(\mathrm{\Lambda }_c^+(2593))M(\mathrm{\Lambda }_c^+)=305.6\pm 0.3`$ MeV and $`h_2^2=0.24_{0.11}^{+0.23}`$, with $`h_2`$ the $`\mathrm{\Lambda }_{c1}\mathrm{\Sigma }_c\pi `$ strong coupling in the chiral Lagrangian. We also comment on the new orbitally excited states recently observed by CLEO.

The charmed baryon system is a convenient testing ground for the ideas and predictions of heavy quark symmetry. This is due to the rich mass spectrum and the relatively narrow widths of the resonances. The properties of these states are the subject of active experimental study at both fixed target experiments (FOCUS, SELEX, E-791) and $`e^+e^{}`$ machines (CLEO, BaBar, Belle). For a recent review of the experimental situation, see Ref. Yelton .

In addition to the usual quantum numbers $`(I,J^P)`$, the charmed baryon states can be labelled also by the spin-parity of the light degrees of freedom $`j_{\mathrm{}}^\pi _{\mathrm{}}`$, which are good quantum numbers in the limit of an infinitely heavy charm quark. This property leads to nontrivial selection rules for the strong couplings of these states to light hadrons IsgWi . These predictions are automatically built into an effective Lagrangian describing the couplings of the heavy baryon states to Goldstone bosons chPT .

The lowest lying charmed baryons are $`L=0`$ states and live in $`\overline{\mathrm{𝟑}}`$ and $`\mathrm{𝟔}`$ representations of flavor SU(3). It is convenient to group them together into superfields defined as in Ref. chPTbar , a vector $`T_i=\frac{1+\mathit{}}{2}(\mathrm{\Xi }_c^0,\mathrm{\Xi }_c^+,\mathrm{\Lambda }_c^+)_i`$ for the $`\overline{\mathrm{𝟑}}`$, and a tensor $`S_\mu ^{ij}=\frac{1}{\sqrt{3}}(\gamma _\mu +v_\mu )\gamma _5B^{ij}+B_\mu ^{ij}`$ for the $`\mathrm{𝟔}`$. These superfields satisfy the constraints from heavy quark symmetry $`\mathit{}T=T`$, $`\mathit{}S_\mu =S_\mu `$ and the condition $`\frac{1+\mathit{}}{2}\gamma ^\mu S_\mu =0`$, which can be used to restrict the form of their Lagrangian interactions sfield . The strong couplings of the lowest lying heavy baryons are described by the effective Lagrangian containing two couplings $`g_{1,2}`$ chPTbar  (we use here the normalization of Ref. PY  for these couplings)

$`_{\mathrm{int}}={\displaystyle \frac{3}{2}}ig_1\epsilon _{\mu \nu \sigma \lambda }(\overline{S}_{ik}^\mu v^\nu A_{ij}^\sigma S_{jk}^\lambda )\sqrt{3}g_2ϵ_{ijk}(\overline{T}^iA_\mu ^{jl}S_{kl}^\mu ),`$ (1)
where $`A_\mu =\frac{i}{2}(\xi ^{}_\mu \xi \xi _\mu \xi ^{})=\frac{1}{f_\pi }_\mu M+\mathrm{}`$ is the usual nonlinear axial current of the Goldstone bosons, defined in terms of $`\xi =\mathrm{exp}(iM/f_\pi )`$ with $`f_\pi =132`$ MeV.

In this paper we focus on the negative parity $`L=1`$ orbitally excited charmed baryons. Combining the quark spins with the $`L=1`$ orbital momentum gives 7 $`\mathrm{\Lambda }`$-type and 7 $`\mathrm{\Sigma }`$-type states without strangeness CIK ; CI  (see Table I). In the constituent quark model, these states fall into two distinct groups, corresponding to the symmetric and antisymmetric irreducible representations of $`S_2`$. The symmetric (antisymmetric) states are denoted in Table I with unprimed (primed) symbols. Quark model estimates for the masses of these states CIK ; CI  suggest that symmetric states are lighter than the antisymmetric ones. Although the permutation symmetry $`S_2`$ is not a true symmetry of QCD beyond the quark model, we will continue to refer to the higher mass charm baryon states as ‘antisymmetric’, as opposed to the lower ‘symmetric’ states. The properties of these states were studied in the quark model in Refs. CIK ; CI ; PY ; JGK  and using large $`N_c`$ methods in CK ; largeNc ; largeNc1 .

The CLEO, ARGUS and E687 Collaborations Lc1  observed two negative parity charm baryons, $`\mathrm{\Lambda }_c^+(2593)`$ and $`\mathrm{\Lambda }_c^+(2625)`$. In accordance with the expectations from the constituent quark model, these states were identified with the $`\mathrm{\Lambda }_{c1}(\frac{1}{2},\frac{3}{2})`$ states in Table I. Their average masses and widths are PDG

$`M(\mathrm{\Lambda }_c^+(2593))M(\mathrm{\Lambda }_c^+)`$ $`=`$ $`308.9\pm 0.6\text{ MeV},\mathrm{\Gamma }(\mathrm{\Lambda }_c^+(2593))=3.6_{1.3}^{+2.0}\text{ MeV}`$ (2)
$`M(\mathrm{\Lambda }_c^+(2625))M(\mathrm{\Lambda }_c^+)`$ $`=`$ $`341.7\pm 0.6\text{ MeV},\mathrm{\Gamma }(\mathrm{\Lambda }_c^+(2625))<1.9\text{ MeV (90\% CL)}`$
| State | $`(I,J)`$ | $`j_{\mathrm{}}^\pi _{\mathrm{}}`$ |
| --- | --- | --- |
| $`\mathrm{\Lambda }_{c1}(\frac{1}{2},\frac{3}{2})`$ | $`(0,\frac{1}{2}),(0,\frac{3}{2})`$ | $`1^{}`$ |
| $`\mathrm{\Sigma }_{c0}(\frac{1}{2})`$ | $`(1,\frac{1}{2})`$ | $`0^{}`$ |
| $`\mathrm{\Sigma }_{c1}(\frac{1}{2},\frac{3}{2})`$ | $`(1,\frac{1}{2}),(1,\frac{3}{2})`$ | $`1^{}`$ |
| $`\mathrm{\Sigma }_{c2}(\frac{3}{2},\frac{5}{2})`$ | $`(1,\frac{3}{2}),(1,\frac{5}{2})`$ | $`2^{}`$ |
| $`\mathrm{\Sigma }_{c1}^{}(\frac{1}{2},\frac{3}{2})`$ | $`(1,\frac{1}{2}),(1,\frac{3}{2})`$ | $`1^{}`$ |
| $`\mathrm{\Lambda }_{c0}^{}(\frac{1}{2})`$ | $`(0,\frac{1}{2})`$ | $`0^{}`$ |
| $`\mathrm{\Lambda }_{c1}^{}(\frac{1}{2},\frac{3}{2})`$ | $`(0,\frac{1}{2}),(0,\frac{3}{2})`$ | $`1^{}`$ |
| $`\mathrm{\Lambda }_{c2}^{}(\frac{3}{2},\frac{5}{2})`$ | $`(0,\frac{3}{2}),(0,\frac{5}{2})`$ | $`2^{}`$ |

> Table I. The quantum numbers of the expected $`p`$-wave strangeless charmed baryons. The corresponding states with strange quarks can be constructed by completing the SU(3) multiplets to which the above states belong.

Motivated by these data, the lowest lying states $`\mathrm{\Lambda }_{c1}(\frac{1}{2},\frac{3}{2})`$ were studied in a chiral Lagrangian approach in Ref. Cho , where their couplings to Goldstone bosons were first derived. These states can be grouped together into a superfield $`R_\mu ^i=\frac{1}{\sqrt{3}}(\gamma _\mu +v_\mu )\gamma _5R^i+R_\mu ^{()i}`$ with $`R_i^{()}=(\mathrm{\Xi }_{c1}^0,\mathrm{\Xi }_{c1}^+,\mathrm{\Lambda }_{c1}^+)_i`$, subject to the same constraints as the superfield $`S_\mu `$.

At leading order in the heavy quark expansion, the pion couplings of these states to the sextet ground state baryons $`S_\mu `$ are given by two terms, corresponding to $`S`$ and $`D`$wave pion emission, respectively

$`_{\mathrm{int}}=h_2ϵ_{ijk}\overline{S}_\mu ^{kl}v_\nu A_{lj}^\nu R_\mu ^i+ih_8ϵ_{ijk}\overline{S}_\mu ^{kl}\left(𝒟^\mu A^\nu +𝒟^\nu A^\mu +{\displaystyle \frac{2}{3}}g^{\mu \nu }(v𝒟)(vA)\right)_{lj}R_\nu ^i+\text{h.c.}`$ (3)
with the covariant derivative $`𝒟_\mu A_\nu =_\mu A_\nu +[V_\mu ,A_\nu ]`$ and $`V_\mu =\frac{1}{2}(\xi ^{}_\mu \xi +\xi _\mu \xi ^{})`$. This formalism was extended to the other $`p`$wave charmed baryons in Table I in Refs. PY ; CF , where prospects were given for their discovery. A total of 6 $`S`$wave and 8 $`D`$-wave couplings are required for a complete description of the strong couplings of the states in Table I.

Knowledge of the pion couplings $`h_2,h_8`$ of the lowest orbital excitations $`\mathrm{\Lambda }_{c1}(\frac{1}{2},\frac{3}{2})`$ will provide information about the other excited baryons, and could thus help guide the search for the missing states. For example, assuming SU(3) symmetry, the widths of the orbitally excited charm baryons containing strange quarks $`\mathrm{\Xi }_{c1}^{}(\frac{1}{2},\frac{3}{2})`$ can be predicted CF ; PY , with results in good agreement with the CLEO data on $`\mathrm{\Xi }_{c1}^{}(\frac{1}{2})`$ CLEOXic1/2  and $`\mathrm{\Xi }_{c1}^{}(\frac{3}{2})`$ CLEOXic3/2 . Furthermore, in the constituent quark model, the couplings of all unprimed states in Table I can be shown to be related to $`h_2,h_8`$ PY ; JGK . Assuming that the masses of these states are known, these relations can be therefore used to predict the decay modes and widths of all these states. Finally, once determined in the charm system, the same couplings would also give the properties of the excited bottom baryons. Clearly, a precise determination of the two couplings $`h_2,h_8`$ is of great interest.

There are a few issues which complicate such a determination, following from the peculiarities of the actual mass spectrum. The states $`\mathrm{\Lambda }_{c1}(\frac{1}{2},\frac{3}{2})`$ are observed through their 3-body decays in the $`\mathrm{\Lambda }_c^+\pi ^+\pi ^{}`$ channel. These are resonant decays, proceeding through intermediate $`\mathrm{\Sigma }_c^{()}\pi `$ states. The masses, and recently the widths of the $`\mathrm{\Sigma }_c`$ baryons have been measured by the FOCUS FOCUS\_sigma  and CLEO CLEO\_sigma  Collaborations. The average results of these measurements are PDG

$`M(\mathrm{\Sigma }_c^{++})M(\mathrm{\Lambda }_c^+)=167.67\pm 0.15\text{ MeV},\mathrm{\Gamma }(\mathrm{\Sigma }_c^{++})=(2.05_{0.38}^{+0.41}\pm 0.38)\text{ MeV}`$ (4)
$`M(\mathrm{\Sigma }_c^+)M(\mathrm{\Lambda }_c^+)=166.4\pm 0.4\text{ MeV},\mathrm{\Gamma }(\mathrm{\Sigma }_c^+)4.6\text{ MeV}(90\%\text{ CL})`$
$`M(\mathrm{\Sigma }_c^0)M(\mathrm{\Lambda }_c^+)=167.32\pm 0.15\text{ MeV},\mathrm{\Gamma }(\mathrm{\Sigma }_c^0)=(1.55_{0.37}^{+0.41}\pm 0.38)\text{ MeV}`$
In the heavy quark limit, the only allowed resonant channels are $`\mathrm{\Lambda }_{c1}(\frac{1}{2})[\mathrm{\Sigma }_c\pi ]_S,[\mathrm{\Sigma }_c^{}\pi ]_D`$, and $`\mathrm{\Lambda }_{c1}(\frac{3}{2})[\mathrm{\Sigma }_c\pi ]_D,[\mathrm{\Sigma }_c^{}\pi ]_{S,D}`$, where the subscript denotes the orbital angular momentum. From (2) and (4) it follows that the dominant $`S`$wave decays of the $`\mathrm{\Lambda }_{c1}(2593)`$ proceed very close to threshold. Furthermore, the available energy in the decay is comparable or less than the width of the decaying state

$`\mathrm{\Lambda }_{c1}(2593)\left[\begin{array}{c}(\mathrm{\Sigma }_c^0(2455)+\pi ^+)\\ (\mathrm{\Sigma }_c^{++}(2455)+\pi ^{})\end{array}\right]\left(\begin{array}{c}2\text{ MeV}\\ 1.7\text{ MeV}\end{array}\right)\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+(2593))`$ (9)
On the other hand, the decay into the $`\mathrm{\Sigma }^+\pi ^0`$ channel takes place $`7.5`$ MeV above threshold, such that it turns out to dominate the width of the $`\mathrm{\Lambda }_{c1}(2593)`$.

The situation with the spin-$`\frac{3}{2}`$ state $`\mathrm{\Lambda }_c(2625)`$ is somewhat different. For this case, the decay is dominated by the $`D`$wave channel $`[\mathrm{\Sigma }_c\pi ]_D`$, which is well above threshold $`(45`$ MeV), while the $`S`$wave accessible modes $`[\mathrm{\Sigma }_c^{}\pi ]_S`$ lie about $`30`$ MeV below threshold and are thus nonresonant.

This suggests that finite width effects are important in the $`\mathrm{\Lambda }_c(2593)`$ decays. The situation is somewhat similar to $`e^+e^{}t\overline{t}`$ production close to threshold, which is mediated by a very broad toponium resonance. The net effect is a distortion of the shape of the invariant mass spectrum in $`\mathrm{\Lambda }_{c1}(2593)\mathrm{\Lambda }_c^+\pi ^+\pi ^{}`$ from a simple Breit-Wigner shape. The resulting line shape depends both on the unknown couplings $`h_{2,8}`$ and on the masses and widths of the intermediate $`\mathrm{\Sigma }_c`$ states. This should be taken into account for the extraction of the mass and width of the $`\mathrm{\Lambda }_{c1}(2593)`$. The purpose of this paper is to present a detailed calculation of these effects.

Consider the amplitude for producing the $`\mathrm{\Lambda }_{c1}`$ resonance, followed by its decay to a 3-body state $`\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi \pi `$, of total momentum $`p_\mu =M_{\mathrm{\Lambda }_c^+}v_\mu +k_\mu `$ and invariant mass $`M(\mathrm{\Lambda }_c^+\pi \pi )=\sqrt{p^2(\mathrm{\Lambda }_c\pi \pi )}`$. This is written in the factorized form

$`A(i\mathrm{\Lambda }_{c1}X\mathrm{\Lambda }_c^+\pi \pi X)={\displaystyle \frac{i}{\mathrm{\Delta }\mathrm{\Delta }_{\mathrm{\Lambda }_{c1}}+i\mathrm{\Gamma }_{\mathrm{\Lambda }_{c1}}(\mathrm{\Delta })/2}}[\overline{U}(\mathrm{\Delta }){\displaystyle \frac{1+\mathit{}}{2}}V(\mathrm{\Delta },X)],`$ (10)
where $`\mathrm{\Delta }=vk=M(\mathrm{\Lambda }_c^+\pi \pi )M(\mathrm{\Lambda }_c^+)`$ is the residual energy of the propagating resonance $`\mathrm{\Lambda }_c(2593)`$ and $`\mathrm{\Delta }_{\mathrm{\Lambda }_{c1}}=M(\mathrm{\Lambda }_{c1})M(\mathrm{\Lambda }_c^+)`$. $`U_\alpha (\mathrm{\Delta })`$ and $`V_\alpha (\mathrm{\Delta },X)`$ are spinor amplitudes parameterizing the decay $`\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi \pi `$ and its production, respectively. $`U_\alpha (\mathrm{\Delta })`$ depends on the momenta and spins of the $`\mathrm{\Lambda }_c\pi \pi `$ state, and is calculable in heavy hadron chiral perturbation theory for values of the residual energy $`\mathrm{\Delta }1`$ GeV. On the other hand, not much is known about the production spinor $`V_\alpha (\mathrm{\Delta },X)`$, which depends on all the details of the production process.

Squaring the amplitude (10), adding the phase space factors and summing over the unobserved states $`X`$, one finds the following expression for the $`\mathrm{\Lambda }_c^+\pi \pi `$ production cross-section as a function of the invariant mass $`\mathrm{\Delta }`$:

$`{\displaystyle \frac{d\sigma (\mathrm{\Delta })}{d\mathrm{\Delta }}}{\displaystyle \frac{1}{(\mathrm{\Delta }\mathrm{\Delta }_{\mathrm{\Lambda }_{c1}})^2+\mathrm{\Gamma }_{\mathrm{\Lambda }_{c1}}(\mathrm{\Delta })^2/4}}\left[\overline{U}(\mathrm{\Delta }){\displaystyle \frac{1+\mathit{}}{2}}\omega (\mathrm{\Delta }){\displaystyle \frac{1+\mathit{}}{2}}U(\mathrm{\Delta })\right]\text{dLips}(\mathrm{\Lambda }_{c1}\mathrm{\Lambda }_c^+\pi \pi )`$ (11)
We have introduced here the density matrix $`\omega _{\alpha \beta }(\mathrm{\Delta })`$ parameterizing the production of a $`\mathrm{\Lambda }_{c1}`$ resonance in the process $`i\mathrm{\Lambda }_{c1}X`$

$`\omega _{\alpha \beta }(\mathrm{\Delta }){\displaystyle \underset{X}{}}{\displaystyle 𝑑\mu (X)V_\alpha (\mathrm{\Delta },X)\overline{V}_\beta (\mathrm{\Delta },X)(2\pi )^4\delta (p_ip_Xp_{\mathrm{\Lambda }_{c1}})}`$ (12)
The matrix $`\omega `$ depends on the resonance momentum $`p_{\mathrm{\Lambda }_{c1}}`$ and details of the experimental setup such as the total beam momentum and polarization. Fortunately, the spin structure of the matrix $`\omega `$ is not required if one sums over the spins and momenta of the final decay products in $`\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi \pi `$. If this is done, the amplitudes in Eq. (11) can be written as

$`{\displaystyle \underset{s_{\mathrm{\Lambda }_c}}{}}{\displaystyle \text{dLips}(\mathrm{\Lambda }_{c1}\mathrm{\Lambda }_c^+\pi ^+\pi ^{})U_\alpha (\mathrm{\Delta })\overline{U}_\beta (\mathrm{\Delta })}=\left({\displaystyle \frac{1+\mathit{}}{2}}\right)_{\alpha \beta }\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi ^+\pi ^{})`$ (13)
Inserting this into (11) one finds that the production cross section as a function of invariant mass takes the factorized form

$`{\displaystyle \frac{d\sigma (\mathrm{\Delta })}{d\mathrm{\Delta }}}\text{Tr }\left[{\displaystyle \frac{1+\mathit{}}{2}}\omega (\mathrm{\Delta })\right]{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi ^+\pi ^{})}{(\mathrm{\Delta }\mathrm{\Delta }_{\mathrm{\Lambda }_{c1}})^2+\mathrm{\Gamma }_{\mathrm{\Lambda }_{c1}}^2(\mathrm{\Delta })/4}}`$ (14)
The dependence on $`\mathrm{\Delta }`$ introduced by the production factor $`\text{Tr }[\frac{1+\mathit{}}{2}\omega (\mathrm{\Delta })]`$ is unknown, and it can be expected to introduce a slow variation with a characteristic scale $`\mathrm{\Lambda }_{QCD}`$. This can be neglected when compared with the much faster variation of the denominator. The width $`\mathrm{\Gamma }(\mathrm{\Delta })`$ in the numerator is equal to the spin-averaged partial width of a $`\mathrm{\Lambda }_{c1}`$ resonance of mass $`\mathrm{\Delta }+M(\mathrm{\Lambda }_c^+)`$ into a specific channel, e.g. $`\mathrm{\Lambda }_c^+\pi ^+\pi ^{}`$, while the width in the denominator $`\mathrm{\Gamma }_{\mathrm{\Lambda }_{c1}}(\mathrm{\Delta })`$ sums over all allowed channels. These decay widths are given explicitly by PY

$`\mathrm{\Gamma }_+\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi ^+\pi ^{})`$ $`=`$ $`{\displaystyle \frac{g_2^2}{16\pi ^3f_\pi ^4}}M_{\mathrm{\Lambda }_c^+}{\displaystyle }dE_1dE_2\{\stackrel{}{p}_2^2|A(E_1,E_2)|^2+\stackrel{}{p}_1^2|B(E_1,E_2)|^2`$ (15)
$`+2\stackrel{}{p}_1\stackrel{}{p}_2\text{Re }[A(E_1,E_2)B^{}(E_1,E_2)]\}`$
where $`E_1,E_2`$ are the pion energies in the rest frame of the $`\mathrm{\Lambda }_{c1}`$ resonance and we have defined

$`A(E_1,E_2)={\displaystyle \frac{h_2E_1}{\mathrm{\Delta }\mathrm{\Delta }_{\mathrm{\Sigma }_c^0}E_1+i\mathrm{\Gamma }_{\mathrm{\Sigma }_c^0}/2}}`$ (16)
$`+h_8\left({\displaystyle \frac{\frac{2}{3}\stackrel{}{p}_1^2}{\mathrm{\Delta }\mathrm{\Delta }_{\mathrm{\Sigma }_c^0}E_1+i\mathrm{\Gamma }_{\mathrm{\Sigma }_c^0}/2}}+{\displaystyle \frac{2\stackrel{}{p}_1\stackrel{}{p}_2}{\mathrm{\Delta }\mathrm{\Delta }_{\mathrm{\Sigma }_c^{++}}E_2+i\mathrm{\Gamma }_{\mathrm{\Sigma }_c^{++}}/2}}\right)`$
$`B(E_1,E_2;\mathrm{\Delta }_{\mathrm{\Sigma }_c^{()0}},\mathrm{\Delta }_{\mathrm{\Sigma }_c^{()++}})=A(E_2,E_1;\mathrm{\Delta }_{\mathrm{\Sigma }_c^{()++}},\mathrm{\Delta }_{\mathrm{\Sigma }_c^{()0}})`$ (17)
The decay rate $`\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi ^0\pi ^0)`$ is given by a similar relation, with an additional factor of $`1/2`$ to account for the identical pions in the final state, and with the replacements $`\mathrm{\Delta }_{\mathrm{\Sigma }_c^{()++}},\mathrm{\Delta }_{\mathrm{\Sigma }_c^{()0}}\mathrm{\Delta }_{\mathrm{\Sigma }_c^{()+}}`$.

In these expressions we work at leading order in the $`1/m_c`$ expansion in matrix elements, but use the exact 3-body phase space. This procedure includes formally subleading contributions in the $`1/m_c`$ expansion, which are however enhanced by kinematics and are required for reproducing the data in other similar situations MLW . We neglect the radiative decay channel $`\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\gamma `$, which is expected to contribute about 20 keV to the total width CK .

After integration over the Dalitz plot, the decay width (15) can be written as

$`\mathrm{\Gamma }_+(\mathrm{\Delta })=g_2^2\left\{h_2^2a_+(\mathrm{\Delta })+h_8^2b_+(\mathrm{\Delta })+2h_2h_8c_+(\mathrm{\Delta })\right\}.`$ (18)
A similar result is obtained for the rate into $`\mathrm{\Lambda }_c^+\pi ^0\pi ^0`$ with coefficients $`a_{00},b_{00},c_{00}`$. The coupling $`g_2`$ appears here both explicitly, and implicitly through the $`\mathrm{\Sigma }_c^{()}`$ widths in the denominators of $`A(E_1,E_2)`$ and $`B(E_1,E_2)`$. These are given by

$`\mathrm{\Gamma }(\mathrm{\Sigma }_c^{()})={\displaystyle \frac{g_2^2}{2\pi f_\pi ^2}}{\displaystyle \frac{M_{\mathrm{\Lambda }_c}}{M_{\mathrm{\Sigma }_c^{()}}}}|\stackrel{}{p}_\pi |^3.`$ (19)
Using the observed masses this gives $`\mathrm{\Gamma }(\mathrm{\Sigma }_c^{++,+,0})=\{6.15,7.06,6.01\}g_2^2`$ MeV, and $`\mathrm{\Gamma }(\mathrm{\Sigma }_c^{++,+,0})=\{47.9,47.4,46.3\}g_2^2`$ MeV. The extracted values for $`g_2`$ from the $`\mathrm{\Sigma }_c`$ and $`\mathrm{\Sigma }_c^{}`$ experimental widths are somewhat different: $`g_2^2_{\mathrm{\Sigma }_c}=0.25\pm 0.17`$, and $`g_2^2_{\mathrm{\Sigma }_c^{}}=0.33\pm 0.15`$, which can be attributed to an $`1/m_c`$ effect. Although the uncertainty in this coupling is rather large, $`g_2^2=0.29\pm 0.23`$, the resulting effect on our predictions (18) is very small, because they are very close to the narrow-width case for the $`\mathrm{\Sigma }_c`$ (see the discussion around Eqs. (20), (21)).

Our main interest here is in the functional dependence of $`a_{+,00}(\mathrm{\Delta })`$, which dominate numerically the rates $`\mathrm{\Gamma }_{+,00}`$. These coefficients are plotted in Fig. 1(a) as functions of $`\mathrm{\Delta }`$; the qualitative features of these curves can be understood without a detailed computation, as follows. The coefficients $`a(\mathrm{\Delta })`$ give the partial widths into the $`[\mathrm{\Sigma }_c\pi ]_S`$ channel, which start at threshold $`\mathrm{\Delta }=2M(\pi ^+)`$, and rise slowly up to the threshold for production of $`[\mathrm{\Sigma }_c^0\pi ^+]_S`$ and $`[\mathrm{\Sigma }_c^{++}\pi ^{}]_S`$ at $`\mathrm{\Delta }=306.9`$ MeV and $`\mathrm{\Delta }=307.2`$ MeV, respectively. Above this threshold, the rate rises much faster, which explains the ‘kink’ seen in Fig. 1(a) in the $`\pi ^+\pi ^{}`$ channel. On the other hand, the threshold in the neutral pion channel lies lower, at $`\mathrm{\Delta }=301.4`$ MeV, corresponding to the opening of the $`[\mathrm{\Sigma }_c^+\pi ^0]_S`$ channel. Since the central value of the $`\mathrm{\Lambda }_{c1}`$ mass lies around $`307`$ MeV, the rapid variation of $`a_+(\mathrm{\Delta })`$ in this region will likely affect the extraction of $`\mathrm{\Delta }_{\mathrm{\Lambda }_{c1}}`$.

It is instructive to compare these results with those obtained in the narrow width approximation, where the mass-dependent partial widths in (15) are approximated with 2-body widths CF

$`\mathrm{\Gamma }_{\mathrm{NW}}(\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi ^+\pi ^{})`$ $`=`$ $`\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+\mathrm{\Sigma }_c^0\pi ^+)+\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+\mathrm{\Sigma }_c^{++}\pi ^{})=a(\pi ^\pm )|\stackrel{}{p}_\pi |`$ (20)
$`\mathrm{\Gamma }_{\mathrm{NW}}(\mathrm{\Lambda }_{c1}^+\mathrm{\Lambda }_c^+\pi ^0\pi ^0)`$ $`=`$ $`\mathrm{\Gamma }(\mathrm{\Lambda }_{c1}^+\mathrm{\Sigma }_c^+\pi ^0)=a(\pi ^0)|\stackrel{}{p}_\pi |`$ (21)
where $`\stackrel{}{p}_\pi `$ is the pion momentum in $`\mathrm{\Lambda }_{c1}\mathrm{\Sigma }_c\pi `$ decays. Neglecting isospin violation in the $`\mathrm{\Sigma }_c`$ masses, the $`a(\pi )`$ parameters are given in the heavy quark limit by

$`a(\pi ^\pm )={\displaystyle \frac{h_2^2}{\pi f_\pi ^2}}{\displaystyle \frac{M_{\mathrm{\Sigma }_c}}{M_{\mathrm{\Lambda }_{c1}}}}E_\pi ^2,a(\pi ^0)={\displaystyle \frac{1}{2}}a(\pi ^\pm ).`$ (22)
In the limit $`g_20`$, the exact result (15) reduces to the narrow width approximation in Eqs. (20) and (21), that is $`\mathrm{\Gamma }\mathrm{\Gamma }_{\mathrm{NW}}`$. As one can see from Fig. 1(a), the narrow width results give a good approximation to the exact widths (computed with $`g_2^2=0.34`$), for $`\mathrm{\Delta }`$ not too close to threshold.

In Fig. 1(b) we show invariant mass distributions $`\mathrm{\Delta }=M(\mathrm{\Lambda }_c^+\pi \pi )M(\mathrm{\Lambda }_c^+)`$ in $`\mathrm{\Lambda }_c^+(2593)`$ decays, in both charged and neutral pions channels. The shape of the invariant mass distribution in the charged pions channel $`\mathrm{\Lambda }_c^+\pi ^+\pi ^{}`$ is distorted towards larger values of $`\mathrm{\Delta }`$ compared to a simple Breit-Wigner curve. In particular, extractions of the $`\mathrm{\Lambda }_c^+(2593)`$ parameters from the charged pions channel alone could overestimate the mass of this resonance by a few MeV, which is larger than the present $`1\sigma `$ uncertainty (2) on this parameter. These effects are not present in the neutral pions channel, for which the shape of the mass spectrum comes closer to a pure Breit-Wigner resonance.

The first observation of the $`\mathrm{\Lambda }_c^+\pi ^0\pi ^0`$ mode has been presented in unpublished CLEO data thesis , where the corresponding invariant mass distribution was used to extract the mass of the $`\mathrm{\Lambda }_c^+(2593)`$. The result is lower than that obtained from the $`\mathrm{\Lambda }_c^+\pi ^+\pi ^{}`$ channel (2), in agreement with our expectations,

$`[M(\mathrm{\Lambda }_c^+(2593))M(\mathrm{\Lambda }_c^+)]_{\mathrm{\Lambda }_c\pi ^0\pi ^0}`$ $`=`$ $`306.3\pm 0.7\text{ MeV}.`$ (23)
Experimental difficulties connected with the low $`\pi ^0`$ detection efficiency could limit the precision of such a determination. We propose therefore that the shape of the $`\mathrm{\Lambda }_c^+\pi ^+\pi ^{}`$ invariant mass spectrum be fit to the distribution (14) with parameters $`(\mathrm{\Delta }_{\mathrm{\Lambda }_{c1}},h_2)`$ (instead of a Breit-Wigner curve with parameters $`(\mathrm{\Delta }_{\mathrm{\Lambda }_{c1}},\mathrm{\Gamma })`$).

In Fig. 2 we show the results of such a fit, performed using the CLEO data presented in thesis  (see Fig. 5.5 in this reference), including detector resolution effects. The parameters of the $`\mathrm{\Lambda }_c(2593)`$ resonance extracted from this fit are <sup>1</sup><sup>1</sup>1The data shown in Fig. 2 was obtained in Ref. thesis  by adding the measured mass difference to a fixed $`\mathrm{\Lambda }_c^+`$ mass of 2286.7 MeV. Thus, for consistency, we subtracted this value from our fitted mass to obtain the result (24).

$`M(\mathrm{\Lambda }_c^+(2593))M(\mathrm{\Lambda }_c^+)=305.6\pm 0.3\text{ MeV},h_2^2=0.24_{0.11}^{+0.23},`$ (24)
and correspond to a resonance mass in reasonably good agreement with (23). A conventional fit of this same data using a Breit-Wigner function, yields a mass difference of around 308 MeV, in agreement with the published measurements Lc1 . Note that the threshold effects effectively lower the resonance mass (24) compared with the previous determinations (2). Our treatment also leads to a reduction in the uncertainties connected with the poorly measured $`\mathrm{\Sigma }_c`$ widths. The result for the coupling $`h_2^2`$ is somewhat lower than previous determinations of this coupling CF  ($`h_2^2=0.30_{0.14}^{+0.21}`$) and PY  ($`h_2^2=0.33_{0.13}^{+0.20}`$).

Finally, we comment on the recent evidence by the CLEO Collaboration CLEO\_new  for new charmed baryon states, lying above the $`\mathrm{\Lambda }_c(2593)`$ and $`\mathrm{\Lambda }_c(2625)`$. The lower signal $`X_1`$ is relatively broad, while the higher peak $`X_2`$ is narrow, with masses and widths

$`M(X_1)M(\mathrm{\Lambda }_c^+)`$ $`=`$ $`480.1\pm 2.4\text{ MeV},\mathrm{\Gamma }(X_1)=20.9\pm 2.6\text{ MeV}`$ (25)
$`M(X_2)M(\mathrm{\Lambda }_c^+)`$ $`=`$ $`595.8\pm 0.8\text{ MeV},\mathrm{\Gamma }(X_2)=4.2\pm 0.7\text{ MeV}`$
The higher narrow peak $`X_2`$ has been identified in CLEO\_new  with the antisymmetric state $`\mathrm{\Lambda }_{c0}^{}(\frac{1}{2})`$, while the lower broad peak $`X_1`$ has been interpreted as the overlap of the peaks corresponding to the two $`\mathrm{\Sigma }_{c2}(\frac{3}{2},\frac{5}{2})`$ states in Table I, which can both decay to $`[\mathrm{\Lambda }_c^+\pi ]_D`$ and $`[\mathrm{\Sigma }_c^{()}\pi ]_D`$. In the heavy quark limit, the $`\mathrm{\Lambda }_{c0}^{}(\frac{1}{2})`$ does not decay to any of the lower lying charmed baryons. The only such mode allowed by isospin and heavy quark symmetry is to $`\mathrm{\Xi }_cK`$, which is however kinematically forbidden. According to the interpretation of CLEO\_new , it can be seen in the $`\mathrm{\Sigma }_c\pi `$ channel through its mixing with the $`\mathrm{\Lambda }_{c1}^{}(\frac{1}{2})`$ at subleading order in $`1/m_c`$. This is consistent with the experimental observation of $`\mathrm{\Sigma }_c\pi `$ resonant substructure (which is accessible in $`\mathrm{\Lambda }_{c1}^{}(\frac{1}{2})[\mathrm{\Sigma }_c\pi ]_S`$), but not of $`\mathrm{\Sigma }_c^{}\pi `$, which can only proceed through $`D`$wave $`\mathrm{\Lambda }_{c1}^{}(\frac{1}{2})[\mathrm{\Sigma }_c^{}\pi ]_D`$ and is therefore expected to be suppressed.

The mass measurement (25) of the $`X_2`$ state shows that it lies above the threshold for $`[ND]`$ ($`\mathrm{\Delta }=524`$ MeV for $`nD^+`$ and $`\mathrm{\Delta }=518`$ MeV for $`pD^0`$). Both the $`\mathrm{\Lambda }_{c0}^{}(\frac{1}{2})`$ and $`\mathrm{\Lambda }_{c1}^{}(\frac{1}{2})`$ can decay to this mode in an orbital $`S`$wave in the heavy quark limit (which would give thus the dominant decay mode in the absence of $`1/m_c`$ effects). Therefore we would like to suggest that one search for the $`X_2`$ state also in the $`ND`$ channel, where it should show up as well. Observing such a signal would definitely rule out alternative interpretations of this state as $`\mathrm{\Sigma }_{c2}(\frac{3}{2},\frac{5}{2})`$ or $`\mathrm{\Lambda }_{c1}^{}(\frac{3}{2}),\mathrm{\Lambda }_{c2}^{}(\frac{3}{2},\frac{5}{2})`$, which can decay to $`ND`$ only in $`D`$wave.

In conclusion, we have discussed in this paper the impact of threshold effects on the determination of the $`\mathrm{\Lambda }_c^+(2593)`$ parameters from its strong decays into $`\mathrm{\Lambda }_c^+\pi \pi `$, and we have presented theory motivated fits of the mass and couplings of this state. Our results suggest that the excitation energy of the $`\mathrm{\Lambda }_c^+(2593)`$ is about 2-3 MeV lower than obtained in previous determinations.

This research was supported by the U.S. NSF Grant PHY-9970781 (A.E.B., A.F.F. and D.P.) and by the DOE grant DE-FG02-97ER41029 (J.M.Y.).


format: text

text:
# Detection of dark matter Skewness in the VIRMOS-DESCART survey: Implications for Ω₀

## 1 Introduction

A direct measurement of the mass distribution in the universe has been a major challenge and focus of modern cosmology. The deflection of light by gravity results in the gravitational lensing effect, which is a direct measure of the strength of gravitational clustering. Recently several groups have been able to measure this weak gravitational lensing effect (Bacon et al., 2002; Refregier et al., 2002; Hoekstra et al., 2002; Van Waerbeke et al., 2002; Jarvis et al., 2002; Brown et al., 2002; Hamana et al., 2002b). Due to the weakness of the effect, all detections have been statistical in nature, primarily in regimes where the signal-to-noise is less than unity.

One of the degeneracies in the measurement of the power spectrum is that between the normalization of the amplitude of fractional fluctuations parameterized by $`\sigma _8`$, the linearly extrapolated standard deviation in spheres of radius $`8h^1`$ Mpc, and the present day density of matter $`\mathrm{\Omega }_0`$. One measures the fluctuations in the projected mass density $`\kappa `$, which can either arise by a large total mass density with small fluctuations, or a smaller mean mass density with larger fractional fluctuations. Typically one constrains the combinations $`\sigma _8\mathrm{\Omega }_0^{0.55}`$(Van Waerbeke et al., 2002). A very similar degeneracy arises in all dynamical measurements of mass, including redshift space distortions (Kaiser, 1987), cluster abundance (Pen, 1998), galaxy peculiar velocities, and pairwise velocities (Davis & Peebles, 1983).

In the standard model of cosmology, fluctuations start off small, symmetric and Gaussian. As fluctuations grow by gravitational instability, this symmetry can no longer be maintained: overdensities can be arbitrarily large, while under dense regions can never have less than zero mass. This leads to a skewness in the distribution of matter fluctuations. The effect has been studied theoretically (Bernardeau et al., 1997) and initial detections have been reported (Bernardeau et al., 2002). In second order perturbation theory, one finds that the skewness scales as the square of the variance, and inversely to density. In terms of the dimensionless surface density $`\kappa `$ and source redshift $`z_s`$, one has

$$S_3\frac{\overline{\kappa }^3}{\overline{\kappa }^2^2}\mathrm{\Omega }_0^{0.8}z_s^{1.35}.$$
(1)
$`S_3`$ does not depend on the power spectrum normalization, but does require knowledge of the source distribution.

The detection of this effect in real data is challenging. The shear field is only sampled at highly irregularly spaced galaxy positions. Just to accurately measure the effective windowed variance of the density field one needs to resort to the two point correlation function (Pen et al., 2002) which can be done cheaply computationally, or through an optimal estimator which requires sophisticated algorithms for the large data sets (Padmanabhan et al., 2002). For low signal-to-noise data, the two point correlation function is the optimal procedure to measure windowed variances and power spectra.

To measure the skewness of a distribution sampled on an irregular pattern requires measuring the three point function. This is itself a computationally complex tasks for the spin 2 weak lensing shear field. Bernardeau et al. (2003) explored a simplified approach using specific geometrical configurations in the shear pattern. They applied their method to the Virmos-Descart data and found the amplitude and the shape of their 3-point shear correlation function over angular scales ranging from one to five arc-minutes follow theoretical expectation of popular cosmological CDM models. Although their detection is strong (4.9 $`\sigma `$), the cosmological interpretation of non-linear features together with the properties of the 2-point shear correlation function is difficult. In particular, one does not know yet how it can be used to break the $`\mathrm{\Omega }_0`$-$`\sigma _8`$ degeneracy, as the skewness of the convergence field can do. Hence, despite recent theoretical investigations of the 3-point shear correlation function (Schneider & Lombardi, 2002; Zaldarriaga & Scoccimarro, 2002; Takada & Jain, 2002), the projection of the three point function onto the skewness of the smoothed convergence is still the easiest and most direct way to constrain cosmological scenario with cosmic shear data.

In this paper we present the first measurement of dark matter skewness from the Virmos-Descart survey111http://terapix.iap.fr/Descart. Despite the strong inhomogeneities in the galaxy distribution produced by the masking, we found an optimal weighting scheme that allows us to compute a robust and reliable skewness estimator. The skewness and variances are computed with the same compensated Gaussian smoothing window, and a full separation into modes has been accomplished. We have performed N-body simulations to calibrate the results.

## 2 Data

The Virmos-Descart data consist in four uncorrelated patches (referred as field F02, F10, F14 and F22 according to their RA position) of about 4 square-degrees each and separated by more than 40 degrees. The fields have been observed with the CFH12k panoramic CCD camera, mounted at the Canada-France-Hawaii Telescope prime focus, over the periods between January 1999 and November 2001.

The observations and data reduction have been described in previous Virmos-Descart cosmic shear papers (Van Waerbeke et al., 2000, 2001, 2002). A detailed discussion on the data quality will be presented elsewhere (McCracken et al in preparation), but most of the data properties relevant for cosmic shear are already presented in Van Waerbeke et al. (2000, 2001, 2002). Shortly, the observations have been done with the I-band filter available on the CFH12k camera with typical exposure time of one hour. The median seeing of the data set is 0.85 arc-second and the limiting depth is $`I_{AB}=24.5`$. All data were processed at the Terapix data center222http://terapix.iap.fr. Although the data observed during that period cover 11.5 sq-degrees, a significant fraction of the field is lost by the masking procedure, so the final area covered by the data discussed in this work is 8.5 sq-degrees.

Catalog construction, galaxy selection, PSF anisotropy correction and shape measurements are presented in Van Waerbeke et al. (2000, 2001, 2002). The final cosmic shear catalog contains 392,055 galaxies with magnitude $`I_{AB}>22`$ with median $`I_{AB}=`$23.6. Several careful checks have demonstrated that systematic residuals are very small. However, Van Waerbeke et al. (2001, 2002) and Pen et al. (2002) have shown that a B-mode signal still remains on scales larger than 10 arc-minutes. Its origin is not yet understood. Below that angular scale, its amplitude is smaller than the E-mode, confirming the gravitational nature of the signal. Van Waerbeke et al. (2001, 2002) and Bernardeau et al. (2002) have used these catalogs to measure the 2-point and 3-points shear correlation functions and constrain cosmological models.

The method can be sumarized as the followingVan Waerbeke et al. (2001): The HDF samples are complete up to the $`I_{AB}=27^{\mathrm{th}}`$ magnitude, and can therefore be used to calibrate the redshift distribution of our VIRMOS sample. In practice, we identify each redshift bin $`n_{m_i}(z)`$ of the HDF’s corresponding to each magnitude bins $`m_i`$ (magnitude bins have a fixed width of $`\mathrm{\Delta }_m=0.5`$). The redshift distribution of the VIRMOS sample is estimated by adding up the $`n_{m_i}(z)`$ with a proper weighting $`w_{m_i}`$. The weight $`w_{m_i}`$ corresponds to the ratio of number counts, per magnitude bin, of the VIRMOS and the HDF catalogues. The completeness of VIRMOS (which occurs 2.5 magnitudes brighter than the HDF) is then properly taken into account.

The resulting histogram for this sample taken from the HDF is shown in Figure 1. In the absence of any spectroscopic survey deeper than $`I_{AB}=22`$ this is the best redshift estimate at the moment. In the future, any accurate cosmological parameter interpretation from the skewness will require precise determination of the source redshift.

## 3 Simulations

### 3.1 N-body

Both the variance and skewness arise from non-linear scales, and simulations are required to calibrate theory (White & Hu, 2000; Jain et al., 2000). The past studies had shown that Gaussian treatments are not accurate for error bar estimates. The sample variance depends on the actual survey geometry, for which simulations are the only tractible approach to quantify these effects.

We ran a series of N-body simulations with different values of $`\mathrm{\Omega }_0`$ to generate convergence maps and make simulated catalogs to calibrate the observational data and estimate errors in the analysis. The power spectra for given parameters were generated using CMBFAST (Seljak & Zaldarriaga, 1996) and these tabulated functions were used to generate initial conditions. The power spectra were normalized to be consistent with the earlier two point analysis from this data set (Van Waerbeke et al., 2002). We ran all of the simulations using a parallel, Particle-mesh N-body code (Dubinski, J., Kim, J., Park, C. 2003, in preparation) at $`1024^3`$ mesh resolution using $`512^3`$ particles on an 8 node quad processor Itanium Beowulf cluster at CITA. Output times were determined by the appropriate tiling of the light cone volume with joined co-moving boxes from $`z6`$ to $`z=0`$. We output periodic surface density maps at $`2048^2`$ resolution along the 3 independent directions of the cube at each output interval. These maps represent the raw output for the run and are used to generate convergence maps in the thin lens and Born approximations by stacking the images with the appropriate weights through the comoving volume contained in the past light cone.

All simulations started at an initial redshift $`z_i=50`$, and ran for 1000 steps in equal expansion factor ratios, with box size $`L=200h^1`$ Mpc comoving. We adopted a Hubble constant $`h=0.7`$ and a scale invariant $`n=1`$ initial power spectrum. A flat cosmological model with $`\mathrm{\Omega }_0+\mathrm{\Lambda }=1`$ was used. Four models were run, with $`\mathrm{\Omega }_0`$ of 0.2, 0.3, 0.4 and 1. The power spectrum normalization $`\sigma _8`$ was chosen as 1.16, 0.90, 0.82 and 0.57 respectively.

We used the simulations in two modes. To study scaling relations and quick analyses, we projected them to convergence maps which were all analyzed with idealized noise on a square domain. We also used the simulations to make mock catalogs which we processed through the full pipeline that was used on the real data. Using the mock catalogs we can quantify the effect of sample variance on $`S_3`$. The errors are expected to be asymmetric: when $`\mathrm{\Omega }_0`$ is large, skewness is small, and the sample variance in skewness is also small. So if one measures a large skewness, one can strongly rule out a high value of $`\mathrm{\Omega }_0`$. But if intrinsic skewness is large, its sample variance is also large. So measuring a small skewness does not rule out a low $`\mathrm{\Omega }_0`$ unless one has a very large field and good statistics.

### 3.2 Simulated Convergence and Shear Maps

Using our N-body simulations we generated two sets of maps. One is a set of idealized $`\kappa `$ convergence maps, from which we can easily measure noise-free statistical quantities. The second is a set of $`\gamma `$ shear maps sampled at the actual 2-D galaxy positions of the survey, which are processed through the same pipeline as the actual observed Virmos-Descart data set. The redshift weights were taken to be the statistical average, so this does not take into account the source-lens-clustering effect (Hamana et al., 2002a; Bernardeau, 1998). Comparing the results from the two procedures allows us to cross-check our analysis, and adds confidence to the interpretation of the complex statistical procedure in terms of simulated viewable dark matter maps.

The convergence $`\kappa `$ is the projection of the matter over-density $`\delta `$ along the line of sight $`\widehat{\theta }`$ weighted by the lensing geometry and source galaxy distribution. It can be expressed as

$$\kappa (\widehat{\theta })=\frac{3}{2}\mathrm{\Omega }_0_0^{\mathrm{}}W(z)\delta (\chi ,\widehat{\theta })𝑑\chi .$$
(2)
$`\chi `$ is the comoving distance in units of $`c/H_0`$. The Hubble constant is parametrized by $`H_0=100h\mathrm{km}/\mathrm{s}/\mathrm{Mpc}`$. The weight function $`W(z)`$ is determined by the source galaxy distribution function $`n(z)`$ and the lensing geometry. For a flat universe,

$$W(z)=(1+z)\chi (z)_z^{\mathrm{}}n(z_s)[1\chi (z)/\chi (z_s)]𝑑z_s.$$
(3)
$`n(z)`$ is normalized such that $`_0^{\mathrm{}}n(z)𝑑z=1`$. For Virmos-Descart, we adopt the distribution shown in figure 1, which is not easily fit by an analytic function.

During each simulation we store 2D projections to the midplane of the dark matter overdensity $`\delta \delta \rho /\overline{\rho }`$ through the 3D box at every light crossing time through the box along each of the x, y and z axes. Our 2D surface density sectional maps are stored on $`2048^2`$ grids. After the simulation, we stack sectional maps separated by the width of the simulation box, randomly choosing the center of each section and randomly rotating and flipping each section. The periodic boundary condition guarantees that there are no discontinuities in any of the maps.

We then add these sections with the weights given by $`W(z)`$ onto a map of constant angular size, which is generally determined by the maximum projection redshift. To minimize the repetition of the same structures in the projection, we alternatively choose the sectional maps of x, y, z directions during the stacking. Using different random seeds for the alignments and rotations, we make $`40`$ maps for each cosmology. Since the galaxy distribution peaks at $`z1`$, the peak contribution of lensing comes from $`z0.5`$ due to the lensing geometry term. Thus a maximum projection redshift $`z2`$ is sufficient for the lensing analysis. We integrate the maps for the $`\mathrm{\Omega }_0=1,0.4,0.3`$ simulations to $`z=2`$ and obtain $`40`$ maps each with angular width $`\theta _\kappa =4.09,3.18`$ and $`3.02`$ degrees, respectively. To make maps for $`\mathrm{\Omega }_0=0.2`$ which are large enough to fit the survey fields, we truncated the integral at to $`z=1.8`$ and obtained maps with angular width $`\theta _\kappa =2.86`$ degree. We show one $`\kappa `$ map from an $`\mathrm{\Omega }_0=0.4`$ simulation in fig. 2. At this resolution, the skewness is quite obvious. The challenge is to recover the skewness from a noisy observation of such a data set.

From the sectional maps, we also produce shear maps. These $`\gamma `$ maps are used to produce mock galaxy catalogs to check the data analysis pipeline and to calibrate the results. For each sectional $`\kappa `$ maps, we computed sectional $`\gamma _1`$ and $`\gamma _2`$ maps using FFT’s. The periodic boundary conditions on each section allow us to circumvent the boundary condition issue. Then by the same projection procedure, we make $`40`$ $`\gamma _1`$ and $`\gamma _2`$ maps with the same fields of view. We sample these simulated shear maps at the actual galaxy position from the Virmos-Descart catalog, which ensures that the mock catalogs have same masks and weights. Since the Virmos-Descart catalog consists of four fields, we can make ten full sets of mock catalogs from each simulation.

With these maps, we simulate the Virmos-Descart survey. Our goal is to measure $`S_3`$ of smoothed $`\kappa `$ fields and its dispersion and distribution in simulated survey catalogs to provide calibrations for the data analysis. In this paper, we fix the smoothing function as the compensated Gaussian function with a filter radius $`\theta _f=5.374`$ arc min. The maps we obtained above are non-periodic after the projection. In order to eliminate the edge effect, we cut $`10\%`$ from the four margins of each smoothed $`\kappa `$ map. For the noise-free maps, the $`S_3`$ parameter is defined as

$$S_3(\theta _f)=\frac{\overline{\kappa }^3}{\overline{\kappa }^2^2}.$$
(4)
We used the $`\kappa `$ maps to fit the dependencies of $`S_3`$ on $`\mathrm{\Omega }_0`$. We parametrized the relation as

$$S_3=\alpha \mathrm{\Omega }_0^\beta .$$
(5)
The best fit values for $`\alpha `$ and $`\beta `$ are

$$\alpha =66.0\pm 2.3,\beta =0.64\pm 0.03,$$
(6)
which are similar to those obtained from second order perturbation theory with different window functions.

## 4 Analysis

### 4.1 Optimality

The optimal estimation of power spectra for Gaussian random fields is well understood. It requires a weighting of the data by the inverse of signal plus noise correlators. This procedure generally costs of order $`O(N^3)`$ where $`N`$ is the number of data points, about one million for our catalog where we have two observables, $`\gamma _1,\gamma _2`$ per galaxy. Clearly the optimal procedure is prohibitively expensive. In principle rapid multi-scale iterative algorithms allow an analysis in $`O(N)`$ (Padmanabhan et al., 2002). The cost prefactor in that algorithm is still large, and we decided to use the two point correlation function instead, whose computational cost is $`O(N\mathrm{log}N)`$ with a small prefactor.

Computing the two point correlation function is in general not an optimal estimator of the power. Being a quadratic estimator, it weights each data bin by only its local noise. Each correlation function bin $`q_i`$ is an inverse noise $`N^1`$ weighted estimator of the data $`\mathrm{\Delta }`$ for a bilinear form $`Q_i`$:

$$q^i=\mathrm{\Delta }^TN^1Q_iN^1\mathrm{\Delta }$$
(7)
while an optimal estimator should have used $`(S+N)^1`$ as its weights instead. In the limit that signal to noise is small, the correlation function is optimal, similarly it is optimal if the signal and noise covariance matrices commute. For weak lensing surveys, on small scales signal to noise is small, while on larger scales the coverage is reasonably uniform, so a correlation function is not a very poor estimator.

To discuss the accuracy of the higher moments of a slightly non-Gaussian random field, consider the one point statistic of a filtered convergence field. The intrinsic distribution of $`\overline{\kappa }`$ is described by $`P_\kappa `$. We observe a noisy realization thereof, $`\overline{\kappa }^o=\overline{\kappa }+\nu `$, with $`\nu `$ denotes a noise random variable. Its variance is $`\sigma _o^2=\sigma _{\overline{\kappa }}^2+\sigma _\nu ^2`$. Typically the variance is noise dominated, in which limit the fractional accuracy of measuring $`\sigma _{\overline{\kappa }}^2`$ when $`\sigma _\nu ^2`$ is known is

$$\frac{\mathrm{\Delta }\sigma _{\overline{\kappa }}^2}{\sigma _{\overline{\kappa }}^2}=\sqrt{\frac{2}{n}}\frac{\sigma _\nu ^2}{\sigma _{\overline{\kappa }}^2}$$
(8)
where $`n`$ is the number of independent samples. On the scales of a few arc minutes where this is best measured, signal to noise is below unity, and an accurate measurement is done statistically where one can reduce the noise due to the large number of arc minute sized patches in the survey.

The skewness is even harder to measure. One expects

$$\frac{\mathrm{\Delta }\overline{\kappa }^3}{\overline{\kappa }^3}=\sqrt{\frac{15}{n}}\left(\frac{\sigma _\nu }{\sigma _{\overline{\kappa }}}\right)^3\frac{1}{S_3\sigma _{\overline{\kappa }}}$$
(9)
The last term is of order unity. Comparing with (8) we notice three losses: 1. The numerical prefactor increases. 2. The error scales as a higher power of signal to noise. 3. The denominator has increased powers of the converge variance, which is itself a small number. These trends suggest that higher moments are increasingly difficult to measure in noisy data.

Optimality of measuring the three point function is a challenging problem. Gaussian fields have zero three point function, so its distribution is intrinsically non-linear and must be measured either through high order perturbation theory or in simulations. By optimality of an estimator we mean one which minimizes the variance. The variance of the three point function is related to the six point function, which is a complex problem.

We thus scale back our ambitions, and concentrate on the skewness in a smoothed density field. In analogy with the two point function, we expect that an inverse noise squared weighting should be a reasonable procedure in this noise dominated measurement.

### 4.2 Two point functions

We first review the efficient computation of two point functions as a warm up exercise for the three point. Computationally, we start with the two shear estimates for each galaxy, and its noise estimate. The weight for each galaxy is the inverse noise variance. The noise variance is the quadratic sum of the intrinsic ellipticity variance $`\sigma _i`$, and the noise variance due to the PSF correction $`\sigma _e`$. The computation of $`\sigma _e`$ is described in Van Waerbeke et al. (2002): the galaxy catalogue is divided into cells of 30 galaxies each in the size-magnitude space. For each cell, the r.m.s. of the ellipticity correction gives an estimate of $`\sigma _e`$. We arbitrarily put a lower bound of 0.28 for the effective noise of any galaxy, as we do not want to give too large a weight to any individual object. Varying this cutoff does not seem to have significant effects on the results. The resulting weight per galaxy is therefore:

$$w_i=\frac{1}{\mathrm{max}(0.28,\sigma _e)^2}.$$
(10)
Galaxies are mapped onto a grid using the cloud-in-cell algorithm (Hockney & Eastwood, 1988). The same is done for the weights. We now have three gridded quantities, which we fast Fourier transform. The autocorrelation of the weights is the inverse Fourier transform of the weight grid times its complex conjugate. Similarly, we obtain the two autocorrelations of the shear components, and two cross correlations, we denote them $`\stackrel{~}{\xi }_{ij}`$. If the dimensions of the grid are a factor of two larger than the spatial extent of the galaxies, the convolution on the grid does not introduce spurious boundary conditions.

The correlation functions are then multiplied by the angular weights, to yield two raw angular averaged correlation functions

$`\stackrel{~}{\xi }_+(r)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle [\xi _{11}(r,\theta )+\xi _{22}(r,\theta )]𝑑\theta }`$
$`\stackrel{~}{\xi }_{}(r)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle [\xi _{11}(r,\theta )\xi _{22}(r,\theta )]\mathrm{cos}(4\theta )}+[\xi _{12}(r,\theta )+\xi _{21}(r,\theta )]\mathrm{sin}(4\theta )d\theta `$ (11)
and an angular averaged weight correlation $`\stackrel{~}{\xi }_w`$. The final correlation is then

$$\xi _\pm =\frac{\stackrel{~}{\xi }_\pm }{\stackrel{~}{\xi }_w}.$$
(12)
This allows a rapid computation of the two point function in eight Fourier transforms, each costing $`2.5N\mathrm{log}_2(N)`$ operations, where $`N=4n_xn_y`$ are the dimensions of the discretized grid, and the factor of 4 comes from having to use a Fourier grid of double the size to deal with the non-periodic boundary conditions. On a PC with the optimized Intel IPP or MKL libraries, the machines can sustain over a gigaflop on a 2 Ghz processor, so using even very fine grids $`n_x10,000`$ takes a matter of seconds.

The variance of the noise is uncorrelated between bins of the two point correlation function. One can see that the covariance of the two point function is a four point function. For a Gaussian field, that can be expanded in terms of two point functions. When the noise is white, a four point function only expands into a non-zero quantity when the four points have two coincident pairs. This is not possible when cross correlating two different correlation lags, and the covariance is thus zero. The error due to noise in each bin of the correlation function is $`4/\stackrel{~}{\xi }_w`$.

In tensor notation, the two point function can be written as a product of functions which depend on the magnitude of separation $`r`$, and angles of the unit separation vector $`\widehat{𝐱}`$

$`\xi _{abcd}(𝐱)`$ $``$ $`\gamma _{ab}(0)\gamma _{cd}(𝐱)`$ (13)
$`=`$ $`{\displaystyle \frac{\xi _+(r)\xi _{}(r)}{2}}(\delta _{ac}\delta _{bd}+\delta _{ad}\delta _{bc}\delta _{ab}\delta _{cd})+4\xi _{}(r)\left(\widehat{x}_a\widehat{x}_b{\displaystyle \frac{\delta _{ab}}{2}}\right)\left(\widehat{x}_c\widehat{x}_d{\displaystyle \frac{\delta _{cd}}{2}}\right)`$
### 4.3 Three point functions

The three point function is computationally more complex. We describe here the algorithm we implemented which computes it in $`O(N^2\mathrm{log}N)`$ operations for $`N`$ galaxies. We first describe the mathematical enumeration of the three point function which is prohibitively expensive. We then describe a faster procedure which computes the same quantity in a shorter time.

At every point we measure two components of the shear. The shear three point function has eight components. To reduce the complexity of the problem, we express the components of the shear in the coordinate system of the line connecting the first two points. We define the first point to be at the origin, and the second point along the x-axis at some distance $`r`$. The third point can be anywhere in the upper half plane. For each galaxy, we first label it point 1. We search rightward, and consider the set of all pairs formed from point 1 and a rightward neighbor, which we successively label as point 2. For each pair 1-2, we consider all galaxies with x coordinates larger than that of galaxy 1, and y coordinates above the line connecting point 1 to point 2. This results in a unique counting of triangles. As for the two point function, we multiply the shear estimates by the weight. We then compute the raw three point functions for the weighted shear, as well as for the weights. The final three point function is the quotient of the two.

At face value this would appear to be an $`O(N^3)`$ operation, which is prohibitively expensive for $`10^5`$ galaxies which we currently have in each field. Analogous to particle-mesh simulations, we map all galaxies onto a two dimensional grid using a chaining mesh. This results in a linked list of galaxies residing in each cell, and costs $`O(N)`$ operations to construct. This allows us to find galaxies within some locus without searching the whole list.

We enumerate all pairs of galaxies. Each pair has a separation and an angle of the line connecting them relative to the x-axis. We grid the list of all pairs onto a four dimensional grid labeled by the position of the first galaxy, and their separations and angles. The grid of separations are chosen logarithmically. The angles run in the range $`(\pi /2,\pi /2)`$, since we only need to count each pair once. For each separation and angle, we convolve this list of pairs with the gridded field of galaxies. We accumulate the three point function by averaging over the angles. In this process, we use only the pair-singlet cross correlation at lags with $`\mathrm{\Delta }x>0`$ and $`y`$ values above the line describing the first pair. This way all triangles are counted exactly once.

In practice, gridding a Virmos-Descart field with 100,000 galaxies on $`512^2`$ cells, using 40 logarithmic separations and 20 angles takes about 15 minutes on the CITA GS320 alpha server which has 32 alpha processors running at 731 Mhz. The procedure requires 15GB of memory to store the four dimensional list of pairs.

The final three point function is defined on a three dimensional grid. In the intermediate pair enumeration we used polar coordinates to save storage and reduce computational cost. Since each pair is tagged by the coordinate of the left galaxy, the three point function has some angular smearing (in our case 9 degrees), so that large relative errors can accrue when galaxies 2 and 3 are nearby to each other but far from galaxy 1. This can be overcome by mapping $`x`$ to $`x`$ and repeating the procedure. For the windowed skewnesses these skinny triangles do not contribute much, and we neglect this effect. In principle one could construct a three point function in $`O(N\mathrm{log}(N))`$ using a tree algorithm. The challenge is that the opening angle has to be chosen small since the shear needs to be rotated by the mean angle.

### 4.4 Skewness

We review the integration of the third moment from the three point function. A smoothed $`\kappa `$ map is defined as

$$\overline{\kappa }(𝐱)=\kappa (𝐱^{})𝒰(|𝐱𝐱^{}|)d^2x^{}$$
(14)
so the third moment can be expressed in terms of an integral over the three point function

$`\overline{\kappa }^3(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{A}}{\displaystyle \kappa (𝐱^{})\kappa (𝐱^{\prime \prime })\kappa (𝐱^{\prime \prime \prime })𝒰(|𝐱𝐱^{}|)𝒰(|𝐱𝐱^{\prime \prime }|)𝒰(|𝐱𝐱^{\prime \prime \prime }|)d^2xd^2x^{}d^2x^{\prime \prime }d^2x^{\prime \prime \prime }}`$ (15)
$`=`$ $`2\pi {\displaystyle \xi _3(r,𝐱)U_3(r,𝐱)r𝑑rd^2x}`$ (16)
where $`A`$ is the area of integration. The three point function $`\xi _3`$ is defined with the first point at the origin, the second along the x axis, and the third anywhere in the plane. The window $`U_3`$ is the overlap integral of three filter functions placed at the same three locations. Due to symmetry, it is sufficient to consider only the upper half of the plane.

The integrands can be evaluated analytically for a compensated Gaussian filter (Crittenden et al., 2002). We choose

$$𝒰=\frac{1}{2\pi r_0^2}\left(1\frac{r^2}{2r_0^2}\right)\mathrm{exp}\left(\frac{r^2}{2r_0^2}+1\right).$$
(17)
This filter has zero area, and is normalized to have a peak amplitude of unity in Fourier space. This has the feature that the filter will only damp modes, and never amplify. Its Fourier transform is $`𝒰(k)=(kr_0)^2\mathrm{exp}(k^2r_0^2/2+1)/2`$. The filter peaks at wave number $`k_0=\sqrt{2}/r_0`$. When $`r_0`$ is measured in radians, $`k`$ is the spherical harmonic number $`l`$. When integrated over a flat spectrum with constant $`l^2C_l`$, the area in Fourier space is $`e^2/8`$. One can use this window as a reasonable broad band Fourier power estimator. One obtains the equivalent flat power by multiplying the variance by $`8/e^2`$. The effective full width at half max width of the filter in Fourier space is a multiplicative factor of 2.34.

The corresponding shear filter is

$$𝒬=𝒰\frac{2}{r^2}_0^rr^{}𝒰(r^{})𝑑r^{}=\frac{r^2}{4\pi r_0^4}\mathrm{exp}\left(\frac{r^2}{2r_0^2}+1\right)$$
(18)
The filtered convergence field is then

$$\overline{\kappa }(𝐱)=\gamma _{ij}(𝐱^{})𝒬(|𝐱𝐱^{}|)\mathrm{\Delta }\widehat{x}_i\mathrm{\Delta }\widehat{x}_jd^2x^{}.$$
(19)
Due to the $`r^2`$ term in front of $`𝒬`$, the hatted objects in the integral are polynomials, and the integrand contains only Gaussian integrals over polynomials, which are exactly solvable. The smoothed variance can be derived in terms of the tensor shear correlation (13)

$$\overline{\kappa }_E^2=\xi _{abcd}(𝐱)W_{abcd}(𝐱)d^2x$$
(20)
for the window function

$`W_{abcd}(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{4(2\pi )^2r_0^4}}{\displaystyle x_a^{}x_b^{}(x_c^{}x_c)(x_d^{}x_d)\mathrm{exp}\left[\frac{x^2+(x^{}x)^2}{2r_0^2}+2\right]d^2x^{}}`$ (21)
$`=`$ $`{\displaystyle \frac{4(\delta _{ac}\delta _{bd}+\delta _{ad}\delta _{bc})r_0^42(\delta _{ac}x_bx_d+\delta _{ad}x_bx_c+\delta _{bc}x_ax_d+\delta _{bd}x_ax_c)r_0^2+x_ax_bx_cx_d}{2^8r_0^4\pi }}`$
$`\times \mathrm{exp}\left({\displaystyle \frac{x^2}{4r_0^2}}+2\right)`$
where we have dropped all trace terms which do not contribute.

For the skewness we can similarly write

$$\overline{\kappa }^3=2\pi \xi _{abcdef}(r,𝐱)W_{abcdef}(r,𝐱)r𝑑rd^2x$$
(22)
where the corresponding window function $`W_{abcdef}`$ is an integral described below.

Since the filter (17) is the Laplacian of a Gaussian, convolving by it is equivalent to first taking the Laplacian of the $`\kappa `$ field, and convolving that by a Gaussian. The Laplacian of the $`\kappa `$ field is also a local second derivative of the shear field, for which a purely local decomposition of the shear field into $`E`$ and $`B`$ mode results independent of survey geometry. We thus expect minimal coupling between $`E`$ and $`B`$ modes from survey geometry.

To construct the analogue for the three point function, it is convenient to define the vectors which point from each vertex to the sum of the other two vertex vectors $`q=x^{}+x^{\prime \prime }`$, $`q^{}=x^{\prime \prime }2x^{}`$ and $`q^{\prime \prime }=x^{}2x^{\prime \prime }`$. They are depicted geometrically in Figure 3.

The non-vanishing terms in the integral are

$`W_{abcdef}`$ $`=`$ $`{\displaystyle \frac{e^3}{(2\pi )^22^33^7r_0^{10}}}[q_aq_bq_c^{}q_d^{}q_e^{\prime \prime }q_f^{\prime \prime }`$ (23)
$`+3\times 4\delta _{ac}q_bq_d^{}q_e^{\prime \prime }q_f^{\prime \prime }+(3\mathrm{cyclic}\mathrm{permutations})`$
$`+9(\delta _{ac}\delta _{be}+\delta _{ae}\delta _{bc})q_d^{}q_f^{\prime \prime }+(12\mathrm{permutations})`$
$`+9\times 2(\delta _{ac}\delta _{bd}+\delta _{ad}\delta _{bc})q_e^{\prime \prime }q_f^{\prime \prime }+(3\mathrm{cyclic}\mathrm{permutations})]`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{x^2+x^{\prime \prime 2}+(x^{\prime \prime }x^{})^2}{6}}\right]`$
To check for systematics, one can compute the moments of the $`B`$ mode skewness $`\overline{\kappa }_B^3`$ by a 45 degree rotation, corresponding to $`(\gamma _1,\gamma _2)(\gamma _2,\gamma _1)`$. Similarly one can get the moments $`\overline{\kappa }_E\overline{\kappa }_B^2`$ and $`\overline{\kappa }_E^2\overline{\kappa }_B`$. When we substitute the three point window (23) into three point integral (16), we place $`𝐱^{}`$ along the x-axis, requiring only a one-dimensional integral, and only integrate $`𝐱^{\prime \prime }`$ on the upper half of the plane, since we only count each triangle once. The sum over the 6 indeces only has 8 independent components, corresponding to the two shear components at each of the three vertices.

In analogy to the two point function, we can compute the expected noise variance in each windowed skewness. The covariance between bins is again zero for white noise, since the only non-trivial six point function consists of three pairs. As for the two point function, we measure the noise from the variance between bins in the three point function, and assume that the noise is linearly proportional to the inverse square root. We verified the procedure by sampling four of the simulated images in the $`\mathrm{\Omega }_0=0.2`$ simulation on a regular $`400^2`$ grid with equal weights, and computing the skewness from the three point function. For these 2.86 degree images this skewness agreed with the image skewness to better than ten percent on the 5.37 arc minute scale, which is consistent with the variation expected from the different weighting of the boundary.

## 5 Results

We show the measured skewness in Figure 4. The variance in the same filter from the two point function is shown in Figure 5. We find that statistical errors from galaxy noise are small, but on all scales we see some $`B`$ mode. The $`B`$ mode is more significant in the variance measurement than in the skewness. This is to be expected, since variances are always positive, while skewness may cancel.

After having obtained these statistics, we were primarily concerned by systematic errors, of which the $`B`$ mode is a potential diagnostic. In our arbitrary binning, the scale at 5.37 arc minutes had the smallest $`B`$ mode in the variance, so we proceeded with all further analyses using this scale. At this scale we find $`\overline{\kappa }^3=1.06\pm 0.06\times 10^6`$. These errors are Gaussian and contain only contributions from noise. Similarly, we find a variance at the same scale of $`\overline{\kappa }^2=6.30\pm 0.62\times 10^5`$. To model the possible effects of the $`B`$ mode, we added the observed $`B`$ power to the error bar, and subtracted the measured $`B`$ power from the signal. Our adopted value for subsequent analysis is $`\overline{\kappa }^2=5.32\pm 0.62\pm 0.98\times 10^5`$, where second error bar is the magnitude of the $`B`$ mode.

The B-mode may or may not contribute to the skewness, which depends on the parity of the source of B modes. The measured value is consistent with zero, so we neglect it in the skewness interpretation. Since variances are always positive, one expects any spurious source of systematics to increase the observed variance. Any corrections due to a observed $`B`$ mode would also correct the observed signal down. For the skewness, the effect could have either sign. We therefore did not correct the $`\overline{\kappa }^3`$ values.

Using the skewness and variances we can construct $`S_3`$. The actual value of the $`B`$ variance is added to to the error bar of the variance. The noise errors are Gaussian since they arise from the mean of a very large number of data points. In the absence of a better model, we also treat the systematic error as Gaussian. The resulting $`S_3`$ is shown in Figure 6. Random variables which are quotients of Gaussians, such as $`S_3`$, don’t generally have well defined moments. The error bars in the plot are obtained by taking one sigma on each of $`\overline{\kappa }^3`$ and $`\overline{\kappa }^2`$, and plotting the observed change in $`S_3`$. For sums of Gaussians this is like summing the standard deviations, which is always an overestimate of the actual standard deviation of the sum. One can consider the error bars in figure 6 to be a conservative estimate. For our scale of 5.37 arc minutes, we used a second error metric by Monte-Carlo sampling of Gaussian noise. Nominally, $`S_3=375`$, and the Monte-Carlo sample finds the median and percentile spreads are $`S_3=376_{125}^{+242}`$ at $`68\%`$ confidence, and $`S_3=376_{194}^{+815}`$ at $`95\%`$. Again, these error bars do not include sample variance, which depends on cosmological models and require N-body simulations to calibrate.

In principle, the fit (5) gives us a value of $`\mathrm{\Omega }_0`$, which would be a very small value, indicating that the data prefers a low density universe with no practical lower bound. The error bar and confidence interval is a tricky task, since it requires knowledge of the area and geometry of the sky patch surveyed. In principle, it should scale as the inverse square root of the number of independent patches as given by Equation (9). At the scales of the filter size of 5.37 arc minutes the survey is quite irregular, and it is not easy to quantify the effective area of sky observed, nor the correlation between patches. So to measure the scatter we use the noise-free mock catalogs from simulations. We process these catalogs through the same three point function simulation pipeline. Even though the simulated catalogs have no noise, we used the same weights as the noisy data. This allows us to separate the error from sample variance from the error from noise.

In figures 7 and 8 we show the results for a $`\mathrm{\Omega }_0=0.4`$ model obtained through the pipeline. The points are the average from ten mock catalogs, and the error bars are the standard deviation across the ten catalogs. The mode separation for the variances is typically better than 1%, while skewness estimates have crosstalk at the 10% level. The plots also include the ensemble averages from the 40 direct idealized $`\overline{\kappa }`$ maps, processed through an independent analysis which shares only the N-body simulation output, which is a useful cross-check of the pipeline. The fit from Peacock & Dodds (1996) is also plotted as the dotted line, which provides a second reference for the resolution of the simulation. We note good agreement for scales larger than 1 arc minute, but below which the simulations are limited by resolution.

While the error bars in figures 4 and 5 only accounted for the noise, the simulation plots show only the sample variance. The sample variance is clearly decreasing at smaller scales, where more independent patches are used. At scales below one arc minute, the simulations are probably limited by resolution.

To estimate the sample variance errors, we used a sample of forty simulated galaxy catalogs. This number is limited by the expense of computations. The sample consisted of our four simulation models with ten catalogs for each model. We scaled the skewness and variance in each model by the inverse of the standard deviation of the same quantity obtained from the simulated maps. This set of normalized samples was used to Monte Carlo sample 10000 instances of $`\overline{\kappa }^3`$ and $`\overline{\kappa }^2`$. For a given trial value of $`\mathrm{\Omega }_0`$, we fixed $`\overline{\kappa }^2`$ to the observed value, and used $`\overline{\kappa }^3`$ from Equation 5. The normalized samples are scaled to these constraints, and we randomly draw from the 40 simulated values. We add noise drawn from a Gaussian with the appropriate variance. This allows us to form a value for $`S_3`$ for each sample. We find the confidence bounds using a uniform Bayesian prior. For $`\mathrm{\Omega }_0=0.25`$ we found 20% of the samples to have $`S_3`$ larger than the observed value, while for $`\mathrm{\Omega }_0=0.5`$ only 10% of the samples hasd $`S_3`$ larger than observed. For $`\mathrm{\Omega }_0=1`$ such large values of $`S_3`$ were never found. We note that our limited simulation sample only provides a crude estimate of the errors. It is probably not meaningful to derive 95% statistics from this small sample whose wings are non-Gaussian. The preliminary interpretation is that a low value of $`\mathrm{\Omega }_0`$ is favored. More detailed systematic numerical studies are in progress. No practical lower bound on $`\mathrm{\Omega }_0`$ is derived, since below $`\mathrm{\Omega }_0=0.2`$ the empirical fit from Equation (5) has not been calibrated.

We compared the sample variance in the simulated maps to that of the mock catalogs. For the scale of interest, we found the scatter in $`\overline{\kappa }^2`$ from the mock catalogs shown in figure 8 consistent with that of a 8.5 sq degree $`\overline{\kappa }`$ map, which is the effective area sampled by our survey. The scatter in $`\overline{\kappa }^3`$ from the mock catalogs was significantly larger than in a 8.5 sq degree contiguous map, possibly closer to what one sees across 3 sq degree maps. This suggests that the inverse noise variance weights are not optimal to estimate skewness. Unfortunately with the limited number of mock catalogs, the error on the sample variance is large, and more simulations are needed. To test the resolution of the simulations, we also compared the results for a simulation with half as many particles in directions on a grid twice as coarse. In the $`\mathrm{\Omega }_0=0.3`$ model, the resulting $`S_3`$ on the $`\overline{\kappa }`$ maps changed by less than one percent. The variance and skewness separately changed by about 10% and 20% respectively, and one would expect the changes to be significantly smaller compared to a higher resolution simulatio. We thus expect that simulation resolution should not be a major source of error on these scales.

We also examined the variance between the four Virmos-Descart fields. At the scale of 5.37 arc minutes, the fields F14, F02, F22 and F10 had $`\overline{\kappa }^3`$ of $`8.7\times 10^8,1.1\times 10^6,5.1\times 10^7,1.1\times 10^6`$ respectively. Their individual values of $`\overline{\kappa }^2`$ were $`3.0\times 10^5,7.0\times 10^5,7.5\times 10^5,7.0\times 10^5`$ without any subtraction of a $`B`$ mode. This scatter is consistent with that seen in the mock catalogs. One could in principle use the observed scatter to model the sample variance, but with only four fields, the error in the scatter is $`\sqrt{2/3}`$ which is very large, so we used the larger sample of simulated mock catalogs.

Several other factors also affect the skewness interpretation. Source-Lens-Clustering (Hamana et al., 2002a; Bernardeau, 1998) for example systematically reduces the observed skewness, which would strengthen the inferred bounds on $`\mathrm{\Omega }_0`$. Using Table 3 from Hamana et al. (2002a), the closest model to our data is probably C1, which for a non-evolving source distribution and the closest comparable window size of 10 arc minutes results in a 16% effect.

The current analysis is limited by sample variance, which will improve dramatically with the CFHT legacy survey. The error on the error estimate is not yet well determined, requiring more simulations to quantify. Errors in the source redshift distribution also effect the results. These effects are expected to be small compared to the sample variance. Work on this issue is in progress, and we expect an order of magnitude improvement in the near future. The precise numerical computation of skewness is also a challenge, as different groups measure slightly different normalizations for the variance (White & Hu, 2000; Jain et al., 2000). For the current data, the differences are not substantial, but need to be tightened with larger and more accurate galaxy catalogs. Currently only one angular scale was used, which is one with the smallest $`B`$ mode contamination in the variance measurement. Even on this optimal scale, the $`B`$ mode will soon become a limiting factor when the survey area is increased.

## 6 Conclusions

We have presented the first direct measurement of dark matter skewness using weak gravitational lensing in the Virmos-Descart survey. We chose the compensated Gaussian with a scale of 5.37 arc minutes for all analyses, which had the smallest $`B`$ mode, consistent with zero. At this scale, the skewness of the dark matter is $`\overline{\kappa }^3=1.06\pm 0.06\times 10^6`$. We measured the variance $`\overline{\kappa }^2=5.32\pm 0.62\pm 0.98\times 10^5`$.

This implies $`S_3=375_{124}^{+342}`$. Calibrating to mock catalogs from N-body simulations, we find $`\mathrm{\Omega }_0<0.5`$ at 90% and $`\mathrm{\Omega }_0<0.25`$ at 80% confidence. The errors are dominated by sample variance. The CFHT legacy survey will increase the sky area by a factor of 20, which should result in a decrease of a factor of 4 in the error in a year’s time. Better constraints could also be obtained by going to smaller scales, where sample variance is smaller. This requires higher resolution simulations, faster three point analysis, and a better understanding of the origin of the $`B`$ mode.

This work was supported in part by the TMR Network “Gravitational Lensing: New Constraints on Cosmology and the Distribution of Dark Matter” of the EC under contract No. ERBFMRX-CT97-0172, the Canada Foundation for Innovation funded computing infrastructure, and NSERC grant 72013704. We thank Francis Bernardeau, Peter Schneider for very interesting and long discussions on skewness and other 3-pt estimators. We thank Uros Seljak for the dark matter projection code, and Mike Jarvis for correction equation (21). TJZ would like to thank CITA for its hospitality during his visit.


text:
# An X-ray Nebula Associated with the Millisecond Pulsar B1957+20

> We have detected an x-ray nebula around the binary millisecond pulsar B1957+20. A narrow tail, corresponding to the shocked pulsar wind, is seen interior to the known H$`\alpha `$ bow shock and proves the long-held assumption that the rotational energy of millisecond pulsars is dissipated through relativistic winds. Unresolved x-ray emission likely represents the shock where the winds of the pulsar and its companion collide. This emission indicates that the efficiency with which relativistic particles are accelerated in the post shock flow is similar to that for young pulsars, despite the shock proximity and much weaker surface magnetic field of this millisecond pulsar.

Millisecond pulsars are old neutron stars (typically $``$ 3 Gyr) that have been spun up to a rapid rotation rate ($`\stackrel{<}{_{}}25`$ ms) by accretion of material from a binary companion (?, ?). After the accretion phase, they appear as radio pulsars with surface magnetic field strengths of $`10^8`$ G, which, combined with their older ages and rapid rotation rates, means that they form a separate population from younger pulsars.

PSR B1957+20 is the second fastest-spinning pulsar known, with a rotation period (?) of 1.6 ms and rotational spin-down luminosity (?) $`\dot{E}=1\times 10^{35}`$ ergs s-1. The pulsar is in a 9.16-hour binary orbit with a low-mass companion star. The wind of the companion star eclipses the radio emission for $``$10% of every orbit. The PSR B1957+20 binary system provides an excellent opportunity to study the wind of a weakly magnetized, recycled neutron star. The wind is ablating, and may eventually evaporate, the low-mass companion star. Ablation and heating of the companion star (?) are believed to be caused by x- or $`\gamma `$-rays generated in an intrabinary shock between the pulsar wind and that of the companion star. Meanwhile, the high space velocity of the pulsar ($`>220`$ km s-1) (?), as it moves through the interstellar medium (ISM), generates sufficient ram-pressure to confine the pulsar wind and results in the formation of a bow shock. Upstream of this bow shock H$`\alpha `$ emission is generated where neutral material is swept up and collisionally excited; this H$`\alpha `$ nebula absorbs 1–10% of the spin-down energy (?). A reverse or termination shock decelerates the pulsar wind and is located between the pulsar and the bow shock.

The intrabinary and termination shocks both provide diagnostics of the pulsar wind. Observations of the interaction of young pulsar winds and the ISM have shown that they power x-ray emitting synchrotron nebulae which are typified by the Crab nebula (?). High-resolution observations of such pulsar-wind nebulae have shown that the composition of the wind, characterized by the ratio of the Poynting flux to particle energy flux $`\sigma `$, is kinetic energy dominated ($`\sigma `$ 0.003). However, for the evolutionarily distinct millisecond pulsar population, the composition, efficiency, and geometry of the pulsar wind remain unknown. ROSAT observations of PSR B1957+20 provided the only previous constraint, and they have been interpreted as indicating that its wind is different from that of the Crab pulsar (?). However, these data lacked the spatial resolution and sensitivity to determine the cause of the emission.

We have undertaken a 43 kilosecond (?) observation of the millisecond pulsar B1957+20 using the Chandra X-ray Observatory. The brightest x-ray source in the field (Fig. 1) is coincident with the pulsar position. A tail of x-ray emission is seen extending from the pulsar to the northeast by at least 16′′ with a position angle opposite to the pulsar’s proper motion direction of 212 (?).

Striking confirmation of the association with the pulsar and the ISM shock comes by comparing this x-ray tail with the H$`\alpha `$ bow-shock nebula (Fig. 2). In bow-shock nebulae, material at the termination shock is swept back by the ram pressure and forms a cylindrical tube aligned with the proper motion direction and interior to the bow shock (?, ?). The morphology, direction and location of the x-ray nebula (Fig. 2) therefore indicate that it corresponds to the termination shock and, in combination with the enclosing H$`\alpha `$ emission, demonstrates the expected double-shock nature of the pulsar’s interaction with the ambient medium.

Young pulsar winds are thought to be relativistic, and correspondingly generate non-thermal emission by synchrotron and/or inverse Compton processes (?, ?). However, nothing is known about millisecond pulsar winds except that the scale of the observed H$`\alpha `$ nebulae are at least consistent with their having relativistic winds (?). A different wind might be expected because of the reduction and possible alteration of the nature of the magnetic field during the accretion phase. We must therefore also consider the possibility that they drive slower winds like that of the Sun, resulting in shock-heated thermal x-rays. If the spin-down luminosity of a millisecond pulsar is carried away as kinetic energy in an outflow then $`\dot{E}=0.5\dot{M}V_w^2`$, where $`\dot{M}`$ is the mass loss rate and $`V_w`$ is the wind speed. Because the pulsar remains and is at least 1 Gyr old we know that $`\dot{M}<10^{17}`$ g s-1, and thus $`V_w>10^9`$ cm s-1. The temperature that we measure ($`kT1`$ keV, where $`k`$ is the Boltzmann constant) by fitting a spectrum expected for shock-heated gas to the energy distribution of the 82 counts recorded in the tail region is inconsistent with such high velocities (?).

The emission must therefore be from some form of non thermal process, either synchrotron or inverse Compton emission, as seen in wind nebulae around young pulsars. In either case, the x-ray emission requires a population of relativistic particles in the pulsar wind. Thus, the detection of a distinct x-ray tail provides direct evidence that millisecond pulsars lose their rotational energy through relativistic winds.

It is highly unlikely that diffuse shock acceleration (?) is the mechanism that generates the emitting particle population in this nebula, because it is difficult to accelerate particles in a relativistic flow through this mechanism (?). As in other pulsar wind nebulae, the acceleration presumably takes place through some other mechanism possibly involving the role of heavy ions (?)

Because the pulsar and the companion are separated by only $`1.5\times 10^{11}`$ cm, the intrabinary shock, formed where the pulsar and companion winds interact, will be located in a strong magnetic field (much stronger than that at the shock in the Crab nebula). This intrabinary shock is therefore a potential source of unresolved synchrotron emission at the location of the pulsar (?) and can be used to determine the wind characteristics at the shock front. The x-ray luminosity in the shock is dependent on both the post shock magnetic field strength and $`\sigma `$ (?). We therefore consider two possibilities for the wind composition of PSR B1957+20: either dominated by kinetic energy ($`\sigma =0.003`$ as seen in the Crab nebula) or magnetically dominated ($`\sigma 1`$). We can describe the process by which the spin-down energy of the pulsar is converted into x-ray emission at the intrabinary shock by (?):

$$f_b\mathrm{\Delta }\epsilon L_\epsilon =f_{rad}f_\gamma f_{geom}\dot{E}$$
(1)
where $`L_\epsilon `$ is the spectral x-ray luminosity in a band of width $`\mathrm{\Delta }\epsilon `$ 1 keV centered on a photon energy $`\epsilon =1`$ keV. The geometric factor $`f_{geom}`$ defines how much of the wind interacts with the companion, $`f_\gamma `$ is the fraction of the intercepted spin-down energy flux that goes into accelerating electrons with Lorentz factor $`\gamma `$ corresponding to $`\epsilon =1`$ keV, and $`f_{rad}`$ is the radiative efficiency of the corresponding synchrotron emission. The fraction of the unresolved x-ray emission that is due to synchrotron emission produced at the intrabinary shock is $`f_b`$, the remainder presumably being emission produced by the pulsar itself.

The flow immediately downstream of the intrabinary shock is expected to undergo Doppler boosting as it passes around the companion (?). Thus, one expects the x-ray emission at orbital phases before and after eclipse to be enhanced by up to a factor 2.2 depending on the flow speed and the degree of absorption and/or scattering by the companion wind (?). In contrast, the x-ray emission at the eclipse (orbital phase 0.25) may be reduced because of obscuration of the shock by the companion star. If we can measure these variations then we can determine $`f_b`$. The lowest and highest count rates in the folded light curve (Fig 3.) are during and immediately after eclipse, respectively, with comparatively low probability of this variation being by chance (see caption of Fig 3.). If this apparent orbital modulation is genuine and corresponds to a modulation of the x-ray flux from the intrabinary shock by a factor of 2.2, then we find that $`f_b0.5`$.

The post shock magnetic field strengths (?) corresponding to our two limiting values of $`\sigma `$ ($`0.003`$ and $`1`$) are listed in Table 1; the corresponding Lorentz factors of the synchrotron emitting relativistic electrons for both cases are $`\gamma =2.4\times 10^5(\epsilon /B)^{1/2}10^5`$ at $`\epsilon =1`$ keV, where $`B`$ is the post shock magnetic field strength. For each case the emitting region is taken to be the radius of the radio-eclipse region (?) $`r_e=5\times 10^{10}`$ cm; for $`\sigma =0.003`$ the flow speed is $`v_{flow}=c/3`$, whereas for $`\sigma 1`$ we expect $`v_{flow}=c`$ (?). The corresponding residence times are given by $`t_{flow}=v_{flow}/r_e`$ (Table 1), and the radiative lifetimes of the emitting electrons are $`t_{rad}=5.1\times 10^8/(\gamma B^2)`$. We can thus compute the radiative efficiency (?) in each case as $`f_{rad}=(1+t_{rad}/t_{flow})^1`$. With $`f_{geom}>0.02`$ (Table 1) and $`f_b0.5`$, we can infer that $`f_\gamma <0.09`$ ($`\sigma =0.003`$) or $`f_\gamma <0.02`$ ($`\sigma 1`$) for Lorentz factor $`\gamma 10^5`$.

We directly compared the properties of the shocked wind of PSR B1957+20 with those derived for the Crab nebula with the same Lorentz factor and find that $`f_\gamma 0.04`$ for electrons with $`\gamma =10^5`$ in the Crab, which, regardless of the assumed value of $`\sigma `$, is similar to that which we derived for PSR B1957+20. The available evidence therefore suggests that despite being subject to a prolonged evolutionary process that has altered its magnetic field by many orders of magnitude and having a shock that occurs much closer to the pulsar, this millisecond pulsar generates a wind for which the efficiency of relativistic particle acceleration in this post shock flow is similar to that seen for the winds of young pulsars. New Chandra data on this and other pulsars (?, ?, ?) are providing the first detailed observational input into studies of relativistic flows, particle acceleration, and magnetohydrodynamic shocks.

* We thank Jon Arons for useful discussions and Heath Jones for providing the H$`\alpha `$ image. This work was supported in part by NASA through a Chandra X-ray Observatory Guest Observer grant, and by the Netherlands Organisation for Scientific Research (NWO).