text:
# Exact time-dependent decoherence factor and its adiabatic classical limit

## I INTRODUCTION

Solvable models in quantum mechanics enable one to investigate quantum measurement problems very conveniently. A good number of investigators have studied these useful models such as Hepp-Coleman model  and Cini model. The exact solvability of these models often provides physicists with a clear understanding of the physical phenomena involved and yield rich physical insights. In these works, the first important step is to obtain the exact solutions of the Schrödinger equation and the time-evolution operator, which can be applied to the calculation of the decoherence factor and study of the wavefunction collapse, etc.. Although the exact solutions and the decoherence of these models have been extensively studied by many authors, the coefficients in the Hamiltonians of all these models are merely time-independent (or partially time-dependent), to the best of our knowledge. In the present paper, we obtain the explicit time-evolution operator and the decoherence factor of the general dynamical models where the Hamiltonians are totally time-dependent.

Time-dependent system is governed by the time-dependent Schrödinger equation. The invariant theory suggested by Lewis and Riesenfeld in 1969 can solve the time-dependent Schrödinger equation. In 1991, Gao et al proposed a generalized invariant theory by introducing basic invariants, which enable one to find the complete set of commuting invariants for some time-dependent multi-dimensional systems. We will analyze the general dynamical model in what follows and then calculate the time-dependent decoherence factor by making use of these invariant theories.

## II EXACT DECOHERENCE FACTOR IN TIME-DEPENDENT DECOHERENCE

The original Cini model for the correlation between the states of the measured system and the measuring instrument-detector is built for a two-level system interacting with the detector. Liu and Sun generalized this Cini model to an $`M`$-level system. In this paper we further consturct an general dynamical model for quantum decoherence between the measuring instrument-detector and the measured system where $`A_+,A_{}`$ and $`A`$ denote the measuring instrument-detector; and $`\omega _i(t)`$ represents the energy parameter of a certain state $`|i`$ of the multi-level measured system and $`\theta _i(t)`$ and $`\varphi _i(t)`$ are coupling coefficients of interaction of the measuring instrument-detector with the measured system. The Hamiltonian which describes the interaction between the state $`|i`$ of the multi-level measured system and the measuring instrument-detector is then given as follows:

$`H_i(t)`$ $`=`$ $`\omega _i(t)\{{\displaystyle \frac{1}{2}}\mathrm{sin}\theta _i(t)\mathrm{exp}[i\varphi _i(t)]A_+`$ (2)
$`+{\displaystyle \frac{1}{2}}\mathrm{sin}\theta _i(t)\mathrm{exp}[i\varphi _i(t)]A_{}+\mathrm{cos}\theta _i(t)A\}`$
with $`A,A_+`$ and $`A`$ satisfying the general commuting relations of a Lie algebra

$$[A_+,A_{}]=nA,[A,A_+]=mA_+,[A,A_{}]=mA_{},$$
(3)
where $`m,m`$ and $`n`$ are structure constants of this Lie algebra. Time evolution of this dynamical model is governed by the Schrödinger equation (in the unit $`\mathrm{}=1`$)

$$i\frac{|\mathrm{\Psi }_i(t)_s}{t}=H_i(t)|\mathrm{\Psi }_i(t)_s.$$
(4)
According to the Lewis-Riesenfeld invariant theory, an operator $`I(t)`$ that agrees with the following invariant equation

$$\frac{I_i(t)}{t}+\frac{1}{i}[I_i(t),H_i(t)]=0$$
(5)
is called an invariant whose eigenvalue is time-independent, i.e.,

$$I_i(t)|\lambda ,i,t=\lambda |\lambda ,i,t_i,\frac{\lambda }{t}=0.$$
(6)
It is seen from Eq. (5) that $`I_i(t)`$ is the linear combination of $`A_+,A_{}`$and $`A`$

and may be generally written

$$I_i(t)=y\{\frac{1}{2}\mathrm{sin}a_i(t)\mathrm{exp}[ib_i(t)]A_++\frac{1}{2}\mathrm{sin}a_i(t)\mathrm{exp}[ib_i(t)]A_{}\}+\mathrm{cos}a_i(t)A,$$
(7)
where the constant $`y`$ will be determined below. Substitution of (7) into Eq. (5) yields

$`y\mathrm{exp}(ib_i)(\dot{a}_i\mathrm{cos}a_ii\dot{b}_i\mathrm{sin}a_i)im\omega _i[\mathrm{exp}(i\varphi _i)\mathrm{cos}a_i\mathrm{sin}\theta _iy\mathrm{exp}(ib_i)\mathrm{sin}a_i\mathrm{cos}\theta _i]`$ $`=`$ $`0,`$ (8)
$`\dot{a}_i+{\displaystyle \frac{ny}{2}}\omega _i\mathrm{sin}\theta _i\mathrm{sin}(b_i\varphi _i)`$ $`=`$ $`0.`$ (9)
where dot denotes the time derivative. The time-dependent parameters $`a_i`$ and $`b_i`$ are determined by these two auxiliary equations.

It is easy to verify that the particular solution $`|\mathrm{\Psi }_i(t)_s`$ of the Schrödinger equation can be expressed in terms of the eigenstate $`|\lambda ,i,t`$ of the invariant $`I_i(t),`$ namely,

$$|\mathrm{\Psi }_i(t)_s=\mathrm{exp}[\frac{1}{i}\phi _i(t)]|\lambda ,i,t$$
(10)
with

$$\phi _i(t)=_0^t\lambda ,i,t^{^{}}\left|[H_i(t^{^{}})i\frac{}{t^{^{}}}]\right|\lambda ,i,t^{^{}}dt^{^{}}.$$
(11)
The physical meanings of $`_0^t\lambda ,i,t^{^{}}\left|H_i(t^{^{}})\right|\lambda ,i,t^{^{}}dt^{^{}}`$ and $`_0^t\lambda ,i,t^{^{}}\left|i\frac{}{t^{^{}}}\right|\lambda ,i,t^{^{}}dt^{^{}}`$ are dynamical and geometric phase, respectively.

Since the expression (10) is merely a formal solution of the Schrödinger equation, in order to get the explicit solutions we make use of the invariant-related unitary transformation formulation which enables one to obtain the complete set of exact solutions of the time-dependent Schrödinger equation (4). In accordance with the invariant-related unitary transformation method, the time-dependent unitary transformation operator is often of the form

$$V_i(t)=\mathrm{exp}[\beta _i(t)A_+\beta _i^{}(t)A_{}]$$
(12)
with $`\beta _i(t)=\frac{a_i(t)}{2}x\mathrm{exp}[ib_i(t)],\beta _i^{}(t)=\frac{a_i(t)}{2}x\mathrm{exp}[ib_i(t)].`$ By making use of the Glauber formula, lengthy calculation yields

$`I_{iV}`$ $`=`$ $`V_i^{}(t)I_i(t)V_i(t)=\{{\displaystyle \frac{y}{2}}\mathrm{exp}(ib_i)\mathrm{sin}a_i\mathrm{cos}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}a_ix]`$ (17)
$`{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{n}}\mathrm{exp}(ib_i)\mathrm{cos}a_i\mathrm{sin}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}a_ix]\}A_+`$
$`+\{{\displaystyle \frac{y}{2}}\mathrm{exp}(ib_i)\mathrm{sin}a_i\mathrm{cos}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}a_ix]`$
$`{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{n}}\mathrm{exp}(ib_i)\mathrm{cos}a_i\mathrm{sin}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}a_ix]\}A_{}`$
$`+\{\mathrm{cos}a_i\mathrm{cos}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}a_ix]+{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}}y\mathrm{sin}a_i\mathrm{sin}[({\displaystyle \frac{mn}{2}})^{\frac{1}{2}}a_ix]\}A.`$
It can be easily seen that when

$$y=\frac{m}{(\frac{mn}{2})^{\frac{1}{2}}},x=\frac{1}{(\frac{mn}{2})^{\frac{1}{2}}},$$
(18)
one may derive that $`I_{iV}=A`$ which is time-independent. Thus the eigenvalue equation of the time-independent invariant $`I_{iV}`$ may be written in the form

$$I_{iV}|\lambda =\lambda |\lambda ,|\lambda =V_i^{}(t)|\lambda ,i,t.$$
(19)
In the meanwhile, by the aid of Baker-Campbell-Hausdorff formula, one can arrive at

$`V_i^{}(t)H_i(t)V_i(t)`$ $`=`$ $`\{\omega _i[{\displaystyle \frac{1}{2}}\mathrm{sin}\theta _i\mathrm{exp}(i\varphi _i){\displaystyle \frac{\sqrt{\frac{mn}{2}}}{n}}\mathrm{exp}(ib_i)\mathrm{cos}\theta _i\mathrm{sin}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)]`$ (24)
$`+{\displaystyle \frac{1}{2}}\omega _i\mathrm{exp}(ib_i)\mathrm{sin}\theta _i\mathrm{cos}(b_i\varphi _i)[\mathrm{cos}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)1]\}A_+`$
$`+\{\omega _i[{\displaystyle \frac{1}{2}}\mathrm{sin}\theta _i\mathrm{exp}(i\varphi _i){\displaystyle \frac{\sqrt{\frac{mn}{2}}}{n}}\mathrm{exp}(ib_i)\mathrm{cos}\theta _i\mathrm{sin}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)]`$
$`+{\displaystyle \frac{1}{2}}\omega _i\mathrm{exp}(ib_i)\mathrm{sin}\theta _i\mathrm{cos}(b_i\varphi _i)[\mathrm{cos}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)1]\}A_{}`$
$`+\omega _i[\mathrm{cos}\theta _i\mathrm{cos}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)+{\displaystyle \frac{\sqrt{\frac{mn}{2}}}{m}}\mathrm{sin}\theta _i\mathrm{cos}(b_i\varphi _i)\mathrm{sin}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)]A,`$
which is related to the dynamical phase; and

$`V_i^{}(t)i{\displaystyle \frac{V_i(t)}{t}}`$ $`=`$ $`[{\displaystyle \frac{i}{2}}\dot{a}_ix\mathrm{exp}(ib_i){\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{\sqrt{\frac{mn}{2}}}}\dot{b}_i\mathrm{exp}(ib_i)\mathrm{sin}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)]A_+`$ (27)
$`+[{\displaystyle \frac{i}{2}}\dot{a}_ix\mathrm{exp}(ib_i){\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{\sqrt{\frac{mn}{2}}}}\dot{b}_i\mathrm{exp}(ib_i)\mathrm{sin}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)]A_{}`$
$`{\displaystyle \frac{\dot{b}_i}{m}}[1\mathrm{cos}(\sqrt{{\displaystyle \frac{mn}{2}}}a_ix)]A,`$
which is related to the geometric phase. It follows from Eq. (24) and Eq. (27) that, under the transformation $`V(t),`$ the Hamiltonian $`H(t)`$ can be changed into

$`H_{iV}(t)`$ $`=`$ $`V_i^{}(t)H_i(t)V_i(t)V_i^{}(t)i{\displaystyle \frac{V_i(t)}{t}}`$ (28)
$`=`$ $`\{\omega _i[\mathrm{cos}a_i\mathrm{cos}\theta _i+{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}}\mathrm{sin}a_i\mathrm{sin}\theta _i\mathrm{cos}(b_i\varphi _i)]`$ (30)
$`+{\displaystyle \frac{\dot{b}_i}{m}}(1\mathrm{cos}a_i)\}A`$
Hence, with the help of Eq. (10) and Eq. (19), the particular solution of the Schrödinger equation is obtained

$$|\mathrm{\Psi }_i(t)_s=\mathrm{exp}[\frac{1}{i}\phi _i(t)]V_i(t)|\lambda ,i$$
(31)
with the phase

$`\phi _i(t)`$ $`=`$ $`{\displaystyle _0^t}\lambda \left|[V_i^{}(t^{^{}})H_i(t^{^{}})V_i(t^{^{}})V_i^{}(t^{^{}})i{\displaystyle \frac{}{t^{^{}}}}V_i(t^{^{}})]\right|\lambda dt^{^{}}`$ (32)
$`=`$ $`\phi _{id}(t)+\phi _{ig}(t)`$ (33)
$`=`$ $`\lambda {\displaystyle _0^t}\{\omega _i[\mathrm{cos}a_i\mathrm{cos}\theta _i+{\displaystyle \frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}}\mathrm{sin}a_i\mathrm{sin}\theta _i\mathrm{cos}(b_i\varphi _i)]`$ (35)
$`+{\displaystyle \frac{\dot{b}_i}{m}}(1\mathrm{cos}a_i)\}\mathrm{d}t^{^{}},`$
where the dynamical phase is $`\phi _{id}(t)=\lambda _0^t\omega _i[\mathrm{cos}a_i\mathrm{cos}\theta _i+\frac{(\frac{mn}{2})^{\frac{1}{2}}}{m}\mathrm{sin}a_i\mathrm{sin}\theta _i\mathrm{cos}(b_i\varphi _i)]dt^{^{}}`$ and the geometric phase is $`\phi _{ig}(t)=\lambda _0^t\frac{\dot{b}_i}{m}(1\mathrm{cos}a_i)dt^{^{}}.`$ It is seen that the former phase is related to the dynamical parameters of the Hamiltonian such as $`\omega _i,\mathrm{cos}\theta _i,\mathrm{sin}\theta _i,`$ etc., whereas the latter is not immediately related to these parameters. If the parameter $`a_i`$ is taken to be time-independent, then we arrive at

$$\phi _{ig}(T)=\lambda _0^T\frac{\dot{b}_i}{m}(1\mathrm{cos}a_i)dt^{^{}}=\frac{\lambda }{m}[2\pi (1\mathrm{cos}a_i)],$$
(36)
where $`2\pi (1\mathrm{cos}a_i)`$ is an expression for the solid angle over the parameter space of the invariant. This fact shows the global or topological meanings of the geometric phase $`\phi _{ig}(t)`$. The expression (31) is a particular exact solution corresponding to $`\lambda `$ and the general solutions of the time-dependent Schrödinger equation are easily obtained by using the linear combinations of all these particular solutions.

Since we have exact solutions of the general time-dependent model, we can obtain the exact expression for the time-dependent decoherence factor that is given

$$F_{i,j}(t)=\lambda \left|V_i^{}(t)V_j(t)\right|\lambda .$$
(37)
Further calculation yields

$$F_{i,j}(t)=\mathrm{exp}[\frac{n\lambda }{2}(\beta _i\beta _j^{}\beta _i^{}\beta _j)]\lambda \left|\mathrm{exp}[(\beta _j\beta _i)A_+(\beta _j^{}\beta _i^{})A_{}]\right|\lambda ,$$
(38)
which is the general expression for the decoherence factor of the time-dependent dynamical model (2). Although the expression (6) is somewhat complicated, it is just the explicit expression that does not contain the chronological product.

## III AN ILLUSTRATIVE EXAMPLE

To show that (6) descends to the result familiar to us in the time-independent (or partially time-dependent) dynamical model of decoherence, we take into consideration the adiabatic classical limit of a special dynamical model. Liu and Sun generalized the original Cini model to an $`M`$-level system. The Hamiltonian of this generalized model is written

$$H=H_S+H_D+H_I,$$
(39)
where $`H_S`$ is the model Hamiltonian of the measured system S with $`M`$ levels and $`H_D`$ is the free Hamiltonian of the two-boson-state detector D. They are generally of the forms

$$H_S=\underset{k=1}{\overset{M}{}}E_k|\mathrm{\Phi }_k\mathrm{\Phi }_k|,H_D=\omega _1a_1^{}a_1+\omega _2a_2^{}a_2$$
(40)
with the creation and annihilation operators $`a_i^{}`$ , $`a_i`$ satisfying the following commuting relations

$$[a_i,a_j^{}]=\delta _{ij},[a_i,a_j]=[a_i^{},a_j^{}]=0.$$
(41)
The interaction Hamiltonian $`H_I`$ is given by

$`H_I`$ $`=`$ $`{\displaystyle \underset{n}{}}|\mathrm{\Phi }_n\mathrm{\Phi }_n|(g_na_1^{}a_2+g_n^{}a_2^{}a_1)`$ (42)
$`=`$ $`{\displaystyle \underset{n}{}}|\mathrm{\Phi }_n\mathrm{\Phi }_n|(g_nJ_++g_n^{}J_{})`$ (43)
with $`J_+=a_1^{}a_2,J_{}=a_2^{}a_1,J_3=\frac{1}{2}(a_1^{}a_1a_2^{}a_2)`$ satisfying the commuting relations $`[J_+,J_{}]=2J_3,[J_3,J_\pm ]=\pm J_\pm .`$

It can be seen from the form of the Hamiltonian that both $`|\mathrm{\Phi }_k\mathrm{\Phi }_k|`$ and $`N=\frac{a_1^{}a_1+a_2^{}a_2}{2}`$ commute with $`H,`$namely, $`[|\mathrm{\Phi }_k\mathrm{\Phi }_k|,H]=[N,H]=0.`$ Hence, a generalized quasialgebra which enables one to obtain the complete set of exact solutions of the Schrödinger equation can be found by working in a sub-Hilbert-space corresponding to the particular eigenvalues of both $`|\mathrm{\Phi }_k\mathrm{\Phi }_k|`$ and $`N`$, and then the Hamiltonian can be rewritten in this sub-Hilbert-space

$$H_{n,k}(t)=E_k+g_kJ_++g_k^{}J_{}+(\omega _1\omega _2)J_3+n(\omega _1+\omega _2)$$
(44)
with $`n`$ being the eigenvalue of $`N`$ and satisfying

$$N|n_1,n_2=n|n_1,n_2,n=\frac{1}{2}(n_1+n_2).$$
(45)
Thus in the sub-Hilbert-space we write the Schrödinger equation in the form

$$H_{n,k}(t)|\mathrm{\Psi }_{n,k}(t)_s=i\frac{}{t}|\mathrm{\Psi }_{n,k}(t)_s,$$
(46)
and $`|\mathrm{\Psi }(t)_s`$ can be obtained from

$$|\mathrm{\Psi }(t)_s=\underset{n}{}\underset{k}{}c_{n,k}|\mathrm{\Psi }_{n,k}(t)_s|\mathrm{\Phi }_k$$
(47)
where $`c_{n,k}`$ is time-independent and determined by the initial conditions.

For the case of this time-dependent Cini model, we set $`A_+=J_+,A_{}=J_{},C=J_z`$ with $`J_\pm =J_1\pm iJ_2`$. In the adiabatic limit, it follows from the auxiliary equations (9) that

$$a_i=\theta _i,b_i=\phi _i;a_j=\theta _j,b_j=\phi _j.$$
(48)
Since

$$\beta _i=\frac{a_i}{2}\mathrm{exp}[ib_i],\beta _j=\frac{a_j}{2}\mathrm{exp}[ib_j],$$
(49)
we let both $`b_i`$ and $`b_j`$ vanish for the convenience, which leads to $`\dot{\theta }_i=0`$ and $`\dot{\theta }_j=0`$ in terms of the auxiliary equations (9), and then $`\beta _i=\frac{a_i}{2},\beta _j=\frac{a_j}{2}.`$ It is therefore easily obtained that

$`F_{i,j}(t)`$ $`=`$ $`\mathrm{exp}[{\displaystyle \frac{\lambda }{2}}(\beta _i\beta _j^{}\beta _i^{}\beta _j)]j,m\left|\mathrm{exp}[(\beta _j\beta _i)J_+(\beta _j^{}\beta _i^{})J_{}]\right|j,m`$ (50)
$`=`$ $`j,m|\mathrm{exp}\{[(\beta _j\beta _i)(\beta _j^{}\beta _i^{})]J_1`$ (52)
$`+i[(\beta _j\beta _i)+(\beta _j^{}\beta _i^{})]J_2\}|j,m`$
$`=`$ $`j,m\left|\mathrm{exp}[i(a_ia_j)J_2]\right|j,m,`$ (53)
where $`j`$ and $`m`$ satisfy

$$J_3|j,m=m|j,m,J^2|j,m=j(j+1)|j,m.$$
(54)
If when $`t=0,`$ then the state of the measuring instrument-detector is $`|j,j,`$ and the decoherence factor is therefore

$`F_{i,j}(t)`$ $`=`$ $`j,j\left|\mathrm{exp}[i(a_ia_j)J_2]\right|j,j`$ (55)
$`=`$ $`[\mathrm{cos}({\displaystyle \frac{a_ia_j}{2}})]^{2j},`$ (56)
which is consistent with that obtained by Sun et al. For the case of classical limit where $`j\mathrm{}`$ and $`\frac{a_ia_j}{2}n\pi (n=0,\pm 1,\pm 2,\mathrm{}),F_{k,l}(t)0`$, which means the wavefunction collapse occurs under the classical limit.

Since we exactly solved the Schrödinger equation governing the time-dependent quantum decoherence, and the result ( i.e., the adiabatic classical limit) obtained here is consistent with what obtained in previous references, we hold that the more general results presented in this paper is useful to treat the time-dependent non-adiabatic quantum decoherence.

## IV CONCLUDING REMARKS

The present paper obtains exact solutions and decoherence factor of the general time-dependent dynamical model for quantum decoherence by working in the sub-Hilbert space corresponding to the eigenvalue of two invariants and by making use of the invariant-related unitary transformation method. The invariant-related unitary transformation formulation is an effective method for treating time-dependent problems. This formulation replaces eigenstates of the time-dependent invariants with those of the time-independent invariants through the unitary transformation. It uses the invariant-related unitary transformation and obtains the explicit expression for the time-evolution operator, instead of the formal solution that is related to the chronological product. In view of what has been discussed above, it can be seen that the invariant theory is appropriate to treat the time-dependent quantum decoherence. Apparently, the results presented in the present paper is easy to generalize to the time-dependent Hepp-Coleman model. Since the geometric phase factor appears in time-dependent systems, it is interesting to investigate the geometric phase in the time-dependent quantum decoherence by using the formulation in this paper.

AcknowledgmentThis project was supported by the National Natural Science Foundation of China under the project No. $`30000034`$.


format: text

text:
# Dark matter distribution in the Coma cluster from galaxy kinematics: breaking the mass-anisotropy degeneracy

## 1 Introduction

The Coma cluster of galaxies (Abell 1656) is one of the most extensively studied in our neighbourhood (see e.g. Biviano 1998 and references therein). Starting with the seminal paper of Kent & Gunn (1982) significant effort went into dynamical modelling of the cluster. In the early studies based on about 300 galaxy velocities only velocity dispersion was modelled and it was most often assumed that the mass follows light and that the galaxies are on isotropic orbits. Merritt (1987) showed that if a larger variety of models is allowed there is a strong degeneracy between the dark matter distribution and velocity anisotropy and many models can be shown to be consistent with the data. Without any prior knowledge on the mass distribution even considering higher velocity moments would probably not be of much help.

Recently, due to theoretical progress mainly by the means of $`N`$-body simulations, our knowledge on possible dark matter distributions within gravitationally bound objects has improved significantly. There seems to be general agreement at least as to the behaviour of dark matter density profiles at large radial distances ($`\varrho r^3`$). Whether the inner dark matter density profile is $`\varrho r^1`$ (as in the so-called universal profile advocated by Navarro, Frenk & White 1997, hereafter NFW) or $`\varrho r^{3/2}`$ (as preferred by Moore et al. 1998; see also Fukushige & Makino 1997) or is flat (as suggested by the observed rotation curves of dwarf and low surface brightness galaxies, e.g. McGaugh & de Blok 1998) is still a matter of debate (a recent analysis by Jimenez, Verde & Oh (2003) of high resolution rotation curves of spiral galaxies shows that 2/3 of the sample can be accounted by NFW profiles, but 2/3 also with a flat core). We therefore consider a generalised profile with different inner slopes and also allow for different dark matter concentrations. We constrain this variety of dark matter density profiles by modelling velocity moments of galaxies. In addition to constraints from the line-of-sight velocity dispersion profile, we incorporate constraints from the fourth velocity moment, the kurtosis.

The NFW profile has been found consistent with the total mass distribution inferred from the galaxy data combined from many clusters in the CNOC1 (van der Marel et al. 2000) and ENACS (Biviano et al. 2003) surveys. Although van der Marel et al. (2000) considered higher velocity moments they did not apply them rigorously to further constrain the mass distribution. Studies based upon X-rays, assuming that the hot X-ray emitting gas is in hydrostatic equilibrium in a spherical potential, usually lead to NFW-like cuspy centres (McLaughlin 1999; Tamura et al. 2000; Sato et al. 2000). The studies based on gravitational lensing focus on the inner shape of the density profile. The slopes agree with the NFW prediction in some studies (e.g. Broadhurst et al. 2000) while one team finds a preference for a flat core (Tyson, Kochanski & dell’Antonio 1998, see also Williams, Navarro & Bartelmann 1999). Note that the Coma cluster is too close for its mass profile to be probed through gravitational lensing.

The amount of galaxy velocity as well as brightness and morphological type measurements for the members of Coma has increased over the last two decades making it possible to analyse separately the samples of elliptical and spiral galaxies. While ellipticals appear to be in dynamical equilibrium justifying the application of Jeans formalism to study their velocity moments, most of the spirals are probably infalling onto the cluster. The modelling of the infall of spirals will be presented in the follow-up study.

The paper is organised as follows. In Section 2 we describe our data. In Section 3 we present our assumptions concerning the matter content of the cluster, i.e. the distributions of galaxies, gas and dark matter. Section 4 outlines our formalism for modelling the velocity moments of elliptical galaxies based on Jeans equations. Section 5 presents its application to constrain the parameters of our model, in particular the distribution of dark matter in the cluster. The discussion follows in Section 6.

## 2 The data

We have searched the NED database for galaxies within $`300^{}`$ of RA=12h59m35.7s, Dec=$`+27^{}57^{}33^{\prime \prime }`$ (J2000) i.e. the position of the elliptical galaxy NGC 4874, well established as the centre of the Coma cluster (Kent & Gunn 1982). The galaxies were required to have velocities between $`3000`$ and $`11000`$ km/s (given the velocity dispersions we will find below, this velocity range extends to $`4\sigma `$). For the calculation of the velocity moments and the subsequent study of kinematics, we remove from the list galaxy pairs and known members of pairs (as given by NED), as we wish to probe the global cluster potential but not its local enhancements. We then obtain a sample of 1068 galaxies shown in the upper panel of Figure 1 as points in the plane of velocity versus distance from the centre of the cluster. To determine the membership of galaxies in the cluster we proceed in a similar fashion as Kent & Gunn (1982). As can be seen from Figure 1, the members of the cluster are well separated from the foreground and background galaxies in velocity space. We have therefore selected the probable members of Coma as lying within the two curves shown in Figure 1, symmetric with respect to $`v=7000`$ km/s, the value close to the mean velocity of the cluster. This procedure leaves us with 967 galaxies. For the determination of the luminosity distribution, we keep the members of galaxy pairs and remove a few galaxies for which no magnitude estimate is available (they may contribute to the calculation of velocity moments however). We then proceed with the membership determination as before.

It is generally believed that only early-type (E-S0) galaxies can be considered in dynamical equilibrium within a cluster in opposition to spirals which are believed to be infalling (e.g. Tully & Shaya 1984; Huchra 1985). We therefore construct separate samples of E-S0 and spiral galaxies. The morphological type of the galaxies has been determined by consulting the NED, SIMBAD and LEDA databases. Among the 967 galaxies belonging to Coma and selected for the analysis of the velocity moments, we find $`355`$ E-S0s, $`163`$ spirals and $`449`$ other galaxies for which the morphological classification is unknown or uncertain. As already noticed by Kent & Gunn (1982), the cluster shows a clear morphological segregation. The distributions of E-S0s and spirals in the $`Rv`$ plane are shown in the middle and lower panel of Figure 1. Comparison of the two plots reveals that while E-S0s cluster and dominate in the central parts of Coma and are not numerous at radial distances larger than $`80^{}`$, spirals are more uniformly distributed and underrepresented in the central region. Similar subsamples are constructed for the analysis of the luminosity distribution. Note that iterating on the mean velocity and establishing new symmetric envelopes for the cluster in the observed phase space has virtually no effect on the cluster membership. We note a possible group of galaxies at $`R<10^{}`$ and $`v<4000`$ km/s that might have been missed, but contributes little to the central internal kinematics.

## 3 The matter content of the cluster

### 3.1 The mass in stars

The mass contributed by the stars in galaxies is estimated as follows. Ideally, one would like to determine the surface luminosity distribution separately for different morphological types of galaxies, since transforming luminosity into mass requires mass-to-light ratios, which are known to vary with morphological type. However, since we do not know the morphological type for roughly half the galaxies in our sample, we would have to assume for these galaxies some mean value of the mass-to-light ratio. Combining then the fits for the three classes of galaxies: E-S0s, spirals and those of unknown type would produce large uncertainties, since the luminosity distribution for spirals turns out to be quite noisy. Therefore, in determining the total stellar mass distribution, we fit the luminosity distribution for all galaxies and then translate it to the mass distribution using a mean mass-to-light ratio.

Magnitudes are transformed into luminosities assuming that all galaxies are at the same distance associated with the mean velocity of the Coma cluster. Adopting an heliocentric velocity of Coma of 6925 km/s (Struble & Rood 1999), and correcting for the velocity of the Sun with respect to the Local Group and for the Local Group infall onto Virgo, we obtain for Coma a Hubble flow velocity of 7093 km/s, which, for a Hubble constant of $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (assumed throughout this paper) gives a distance of 101.3 Mpc and a distance modulus of 35.03 (neglecting the peculiar velocity of the Coma cluster).

As mentioned in the previous section, our sample of galaxies for the luminosity analysis is somewhat different from the one used for the calculation of velocity moments, as we no longer exclude the members of pairs. Now, we have 985 galaxies that belong to Coma, among which 366 E-S0s and 167 spirals. The surface luminosity profile of all galaxies is then determined by placing 60 galaxies per radial bin. The resulting distribution is shown as filled symbols in Figure 2. The open symbols in the Figure show a similar result for just E-S0 galaxies but this time with 30 galaxies per radial bin. This second distribution will be needed in the modelling of the velocity moments as the distribution of our tracer population.

The data shown in Figure 2 do not indicate any presence of a core in the surface luminosity distributions, hence we fit them with projections of cuspy profiles. For both samples, of all galaxies as well as early-type ones, the distributions have a changing slope, so we fit them with a projection of the NFW profile (Bartelmann 1996; see also section 2.5 of Łokas & Mamon 2001). The fit is done by $`\chi ^2`$ minimisation taking all points with equal weights (although we do not know the errors of the magnitude measurements we expect them to be similar for galaxies in each bin). The 3D luminosity density will then be

$$l(r)=\frac{l_{}}{(r/r_\mathrm{S})(1+r/r_\mathrm{S})^2},$$
(1)
where the normalisation constant $`l_{}`$ and the scale radius $`r_\mathrm{S}`$ are the two fitting parameters. For the sample of all galaxies we obtain $`l_{}=9.05\times 10^7L_{\mathrm{}}/`$arcmin3 and $`r_\mathrm{S}=14.^{}4`$. For the E-S0s we have $`l_{}=3.55\times 10^8L_{\mathrm{}}/`$arcmin3 and $`r_\mathrm{S}=7.^{}05`$ so this distribution turns out to be more concentrated and therefore steeper. It might be interesting to note that the *number* density distribution of the galaxies in our sample is less steep in the centre than the *luminosity* density distribution and can be approximated by a profile with a core. The steeper distribution of luminosity is probably mainly due to the presence of bright ellipticals in the centre of the cluster (see Figure 1).

Integrating the luminosity distribution (1) and multiplying by the appropriate mass-to-light ratio in the blue band, $`\mathrm{{\rm Y}}`$, we obtain the mass distribution associated with the stars in galaxies

$$M_\mathrm{G}(r)=4\pi L_{}\mathrm{{\rm Y}}r_\mathrm{S}^3\left(\mathrm{ln}\frac{r+r_\mathrm{S}}{r_\mathrm{S}}\frac{r}{r+r_\mathrm{S}}\right).$$
(2)
We adopt the mass-to-light ratio in blue band for E-S0s of $`\mathrm{{\rm Y}}_\mathrm{E}=8M_{\mathrm{}}/L\mathrm{}`$ and for spirals $`\mathrm{{\rm Y}}_\mathrm{S}=3M_{\mathrm{}}/L_{\mathrm{}}`$. For the mass distribution of all galaxies we take the weighted mean of $`\mathrm{{\rm Y}}`$. Since, for galaxies with known morphological type, the E-S0s and spirals appear in the proportions of $`2.2:\mathrm{\hspace{0.17em}1}`$, then assuming that the same proportion holds for the whole population, the mean mass-to-light ratio amounts to $`\mathrm{{\rm Y}}_\mathrm{G}=(2.2\mathrm{{\rm Y}}_\mathrm{E}+\mathrm{{\rm Y}}_\mathrm{S})/3.2=6.43M_{\mathrm{}}/L_{\mathrm{}}`$.

### 3.2 The gas distribution

To approximate the contribution of the gas to the mass in the cluster we make use of the X-ray surface brightness distribution which can be well approximated by the following formula

$$S(r)=S_0\left[1+\left(\frac{r}{r_\mathrm{c}}\right)^2\right]^{3b+1/2}.$$
(3)
Briel, Henry & Böhringer (1992) analysed the ROSAT observations of Coma out to $`100^{}`$ in the energy band of 0.5 to 2.4 keV and obtained the following best-fitting parameters: $`S_0=4.6\times 10^{13}`$ erg cm-2 s-1 arcmin-2, $`b=0.75`$ and the core radius $`r_\mathrm{c}=10.^{}5`$. Assuming that the radiation is produced by bremsstrahlung this distribution is obtained by projecting the 3D emissivity integrated over the energy band and assuming that the gas is approximately isothermal. Since the emissivity is proportional to the square of the number density of electrons in the gas $`n`$, we obtain

$$n(r)=n_0\left[1+\left(\frac{r}{r_\mathrm{c}}\right)^2\right]^{3b/2},$$
(4)
where the central electron density (with $`h=0.7`$) is $`n_0=3.42\times 10^3`$ cm-3 (Briel et al. 1992; Henry & Henriksen 1986). Integrating equation (4) and multiplying by the mass per electron we obtain the mass distribution associated with the gas

$$M_\mathrm{g}(r)=\frac{4}{3}\pi n_0(m_\mathrm{e}+\gamma m_\mathrm{p})r_2^3F_1(\frac{3}{2},\frac{3b}{2};\frac{5}{2};\frac{r^2}{r_\mathrm{c}^2}),$$
(5)
where $`{}_{2}{}^{}F_{1}^{}`$ is the hypergeometric function and $`\gamma m_\mathrm{p}`$ with $`\gamma =1.136`$ is the mean mass of the positively charged ion in the gas per unit charge ($`m_\mathrm{p}`$ being the proton mass), assuming the gas has the primordial composition with helium abundance of $`Y_{\mathrm{He}}=0.230.24`$.

The formula (3) for the surface brightness is the same as the one appearing in the so-called $`\beta `$-model (we have replaced $`\beta `$ by $`b`$ to avoid confusion with the anisotropy parameter of the next section), however we do not accept all the assumptions of the model here concerning e.g. the relation between the gas and galaxy distributions. The only assumption going into the derivation of the gas mass distribution in equation (5), besides spherical symmetry, is that the gas is isothermal, which is justified by recent observations of Coma with XMM-Newton by Arnaud et al. (2001), who find very little temperature variation at least in the central region of radius $`20^{}`$.

### 3.3 Dark matter

We study different possible density distributions of dark matter described by the following general formula

$$\rho (r)=\frac{\rho _{\mathrm{char}}}{(r/r_\mathrm{s})^\alpha (1+r/r_\mathrm{s})^{3\alpha }},$$
(6)
where $`\rho _{\mathrm{char}}`$ is a constant characteristic density, while $`r_\mathrm{s}`$ is the scale radius of the dark matter (in general different from that of the luminous distribution, $`r_\mathrm{S}`$). The profiles differ by the inner slope $`r^\alpha `$ but have a common outer limiting behaviour of $`r^3`$. We will consider $`\alpha `$ values limited by $`0\alpha 3/2`$, which covers a wide range of possible inner profiles: from very steep to core-like. The cuspy profiles of $`\alpha >0`$ are motivated by the results of cosmological $`N`$-body simulations. The profile with $`\alpha =1`$ corresponds to the so-called universal profile proposed by NFW as a fit to the profiles of simulated haloes, while the profile with $`\alpha =3/2`$ is identical to the one following from higher resolution simulations of Moore et al. (1998). The core profile with $`\alpha =0`$ is favoured by some observations of galaxies and clusters and is very similar (but not identical) to the profile proposed by Burkert (1995).

The scale radius $`r_\mathrm{s}`$ introduced in equation (6) marks the distance from the centre of the object where the profile has a slope equal to the average of the inner and outer slope: $`r^{(3+\alpha )/2}`$. The other parameter that controls the shape of the profile is the concentration

$$c=\frac{r_v}{r_\mathrm{s}},$$
(7)
where $`r_v`$ is the virial radius, i.e. the distance from the centre of the halo within which the mean density is $`\mathrm{\Delta }_\mathrm{c}`$ times the present critical density, $`\rho _{\mathrm{crit},0}`$. Although in most cosmological $`N`$-body simulations $`\mathrm{\Delta }_\mathrm{c}=200`$ is assumed and kept constant, the value, following from the spherical collapse model, depends on the underlying cosmology, e.g. $`\mathrm{\Delta }_\mathrm{c}178`$ is valid for the Einstein-de Sitter model, while, for the currently most popular $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_M=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, we have $`\mathrm{\Delta }_\mathrm{c}102`$ (Eke, Cole & Frenk 1996; Łokas & Hoffman 2001). We will keep $`\mathrm{\Delta }_\mathrm{c}=102`$ in the following.

The concentration of simulated dark matter haloes has been observed to depend on the virial mass. Jing & Suto (2000) tested the relation $`c(M_v)`$ for the masses of the order of normal galaxies and clusters in the case of density profiles with $`\alpha =1`$ and $`\alpha =3/2`$ and found concentrations slowly decreasing with mass (thus confirming the original observation of NFW) and lower for $`\alpha =3/2`$ than for $`\alpha =1`$. The only study using the value of $`\mathrm{\Delta }_\mathrm{c}`$ appropriate for a given cosmological model is the one of Bullock et al. (2001) who found the profiles of presently forming haloes to be well fitted by the NFW formula with concentrations depending on mass in the $`\mathrm{\Lambda }`$CDM model with the above parameters approximately as

$$c(M_v)=5.95\left(\frac{M_v}{10^{15}h^1M_{\mathrm{}}}\right)^{0.122}.$$
(8)
In the following, we will treat the concentration as a free parameter using this relation to guide us as to order of magnitude of $`c`$.

We normalise the density profile (6) so that the mass within $`r_v`$ is equal to the so-called virial mass

$$M_v=\frac{4}{3}\pi r_v^3\mathrm{\Delta }_\mathrm{c}\rho _{\mathrm{crit},0}.$$
(9)
The characteristic density of equation (6) then becomes

$$\rho _{\mathrm{char}}=\frac{(3\alpha )\mathrm{\Delta }_\mathrm{c}\rho _{\mathrm{crit},0}c^\alpha }{3F_\alpha (c)},$$
(10)
where $`F_\alpha (c)`$ is given by the hypergeometric function

$`F_\alpha (x)`$ $`=`$ $`{}_{2}{}^{}F_{1}^{}(3\alpha ,3\alpha ;4\alpha ;x)`$ (11)
$`=`$ $`\{\begin{array}{cc}{\displaystyle \frac{3}{x^3}}\left[\mathrm{ln}(1+x){\displaystyle \frac{x}{2}}{\displaystyle \frac{(2+3x)}{(1+x)^2}}\right]\hfill & (\alpha =0)\hfill \\ {\displaystyle \frac{5}{x^{5/2}}}\left[\mathrm{sinh}^1\sqrt{x}{\displaystyle \frac{\sqrt{x}}{3}}{\displaystyle \frac{3+4x}{(1+x)^{3/2}}}\right]\hfill & (\alpha =1/2)\hfill \\ {\displaystyle \frac{2}{x^2}}\left[\mathrm{ln}(1+x){\displaystyle \frac{x}{1+x}}\right]\hfill & (\alpha =1)\hfill \\ {\displaystyle \frac{3}{x^{3/2}}}\left[\mathrm{sinh}^1\sqrt{x}\sqrt{{\displaystyle \frac{x}{1+x}}}\right]\hfill & (\alpha =3/2)\hfill \end{array}`$ (16)
The dark mass distribution following from (6), (9) and (10) is

$$M_\mathrm{D}(s)=M_vs^{3\alpha }\frac{F_\alpha (cs)}{F_\alpha (c)},$$
(17)
where we introduced $`s=r/r_v`$.

## 4 Modelling of the velocity moments

The purpose of this work is to constrain the distribution of dark matter in the Coma cluster by studying the velocity moments of the population of elliptical galaxies in the cluster which we believe to be in equilibrium (thus neglecting radial streaming motions). We defer to a later paper the kinematical study of spiral galaxies in the context of infall. We will also assume that the system is spherically symmetric and that there are no net streaming motions (e.g. no rotation) so that the odd velocity moments vanish.

At second order, the two distinct moments are $`\overline{v_r^2}`$ and $`\overline{v_\theta ^2}=\overline{v_\varphi ^2}`$ which we will denote hereafter by $`\sigma _r^2`$ and $`\sigma _\theta ^2`$ respectively. They can be calculated from the lowest order Jeans equation (e.g. Binney & Mamon 1982)

$$\frac{\mathrm{d}}{\mathrm{d}r}(\nu \sigma _r^2)+\frac{2\beta }{r}\nu \sigma _r^2=\nu \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r},$$
(18)
where $`\nu `$ is the 3D density distribution of the tracer population and $`\mathrm{\Phi }`$ is the gravitational potential. We will solve equation (18) assuming the anisotropy parameter

$$\beta =1\frac{\sigma _\theta ^2(r)}{\sigma _r^2(r)}$$
(19)
to be constant with $`\mathrm{}<\beta 1`$. This model covers all interesting possibilities from radial orbits ($`\beta =1`$) to isotropy ($`\beta =0`$) and circular orbits ($`\beta \mathrm{}`$).

The solution of the lowest order Jeans equation with the boundary condition $`\sigma _r0`$ at $`r\mathrm{}`$ for $`\beta `$=const is

$$\nu \sigma _r^2(\beta =\mathrm{const})=r^{2\beta }_r^{\mathrm{}}r^{2\beta }\nu \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}dr.$$
(20)
However, the measurable quantity is the line-of-sight velocity dispersion obtained from the 3D velocity dispersion by integrating along the line of sight (Binney & Mamon 1982)

$$\sigma _{\mathrm{los}}^2(R)=\frac{2}{I(R)}_R^{\mathrm{}}\left(1\beta \frac{R^2}{r^2}\right)\frac{\nu \sigma _r^2r}{\sqrt{r^2R^2}}dr,$$
(21)
where $`I(R)`$ is the surface distribution of the tracer and $`R`$ is the projected radius. Introducing equation (20) into equation (21) and inverting the order of integration, we obtain

$`\sigma _{\mathrm{los}}^2(R)`$ $`=`$ $`{\displaystyle \frac{2G}{I(R)}}{\displaystyle _R^{\mathrm{}}}dx\nu (x)M(x)x^{2\beta 2}`$
$`\times {\displaystyle _R^x}\mathrm{d}y(1\beta {\displaystyle \frac{R^2}{y^2}}){\displaystyle \frac{y^{2\beta +1}}{\sqrt{y^2R^2}}},`$
where $`M(x)`$ is the mass distribution and we used new variables $`x`$ and $`y`$ instead of $`r`$ to avoid confusion. The calculation of $`\sigma _{\mathrm{los}}`$ can then be reduced to one-dimensional numerical integration of a formula involving special functions for arbitrary $`\beta =\mathrm{const}`$.

It has been established that by studying $`\sigma _{\mathrm{los}}(R)`$ alone we cannot uniquely determine the properties of a stellar system. In fact, systems with different densities and velocity anisotropies can produce identical $`\sigma _{\mathrm{los}}(R)`$ profiles (see e.g. Merrifield & Kent 1990; Merritt 1987). It is therefore interesting to consider higher-order moments of the velocity distribution. For the fourth-order moments, the three distinct components $`\overline{v_r^4}`$, $`\overline{v_\theta ^4}=\overline{v_\varphi ^4}`$ and $`\overline{v_r^2v_\theta ^2}=\overline{v_r^2v_\varphi ^2}`$ are related by two higher order Jeans equations (Merrifield & Kent 1990).

In order to solve these equations we need additional information about the distribution function. We will restrict ourselves here to functions which can be constructed from the energy-dependent distribution function by multiplying it by a function of angular momentum $`f(E,L)=f_0(E)L^{2\beta }`$ with $`\beta =\mathrm{const}`$. The solution of the Jeans equation for the fourth-order moment (see Łokas 2002)

$$\frac{\mathrm{d}}{\mathrm{d}r}(\nu \overline{v_r^4})+\frac{2\beta }{r}\nu \overline{v_r^4}+3\nu \sigma _r^2\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}=0,$$
(23)
is

$$\nu \overline{v_r^4}(\beta =\mathrm{const})=3r^{2\beta }_r^{\mathrm{}}r^{2\beta }\nu \sigma _r^2(r)\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}dr.$$
(24)
By projection, we obtain the line-of-sight fourth moment

$$\overline{v_{\mathrm{los}}^4}(R)=\frac{2}{I(R)}_R^{\mathrm{}}\frac{\nu \overline{v_r^4}r}{\sqrt{r^2R^2}}g(r,R,\beta )dr,$$
(25)
where

$$g(r,R,\beta )=12\beta \frac{R^2}{r^2}+\frac{\beta (1+\beta )}{2}\frac{R^4}{r^4}.$$
(26)
Introducing equations (20) and (24) into (25) and inverting the order of integration, we obtain

$`\overline{v_{\mathrm{los}}^4}(R)`$ $`=`$ $`{\displaystyle \frac{6G^2}{I(R)}}{\displaystyle _R^{\mathrm{}}}{\displaystyle \frac{r^{2\beta +1}}{\sqrt{r^2R^2}}}g(r,R,\beta )dr`$ (27)
$`\times {\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{\nu (q)M(q)}{q^{22\beta }}}\mathrm{d}q{\displaystyle _r^q}{\displaystyle \frac{M(p)}{p^2}}\mathrm{d}p`$
and the calculation of $`\overline{v_{\mathrm{los}}^4}(R)`$ can be reduced to a (rather long) double integral. A useful way to express the fourth projected moment is to scale it with $`\sigma _{\mathrm{los}}^4`$ in order to obtain the projected kurtosis

$$\kappa _{\mathrm{los}}(R)=\frac{\overline{v_{\mathrm{los}}^4}(R)}{\sigma _{\mathrm{los}}^4(R)},$$
(28)
whose value is 3 for a Gaussian distribution.

## 5 Results

We use our sample of galaxy velocities to calculate the velocity moments of the E-S0 galaxies. The 355 available E-S0 galaxies were divided into bins of 39 objects. Figure 3 shows the line-of-sight velocity dispersion (upper panel) and kurtosis (lower panel) together with their sampling errors calculated using estimators based on Monte Carlo simulations (see Appendix A). We find the sampling distribution of $`\sigma _{\mathrm{los}}`$ to be close to normal. In the case of kurtosis, shown in the lower panel of Figure 3 the values are actually calculated and their errors propagated from the quantity $`(\mathrm{log}\kappa _{\mathrm{los}})^{1/10}`$ which we find to be normally distributed (see Appendix A). Normal sampling distributions of the estimators of both moments and very weak correlations between them justify the use of standard fitting procedures of these quantities based on $`\chi ^2`$ minimisation.

Our purpose here is to reproduce the observed velocity moments using models described by equations (21) and (27), with the mass distribution given by the sum of the three contributions (2), (5) and (17) discussed in Section 3:

$$M(r)=M_\mathrm{G}+M_\mathrm{g}+M_\mathrm{D}.$$
(29)
The density profile $`\nu (r)`$ of the tracer population of early-type galaxies is given by equation (1) and the surface brightness $`I(R)`$ by its projection. While studying velocity dispersion is useful to constrain the mass, the kurtosis is mostly sensitive to the velocity anisotropy.

To give a feeling of its behaviour, we show in Figure 3 the predictions of equation (27) (and the corresponding ones of equation (21)) for the three cases of radial ($`\beta =1`$), isotropic ($`\beta =0`$) and circular ($`\beta =\mathrm{}`$) orbits assuming dark matter distribution given by an NFW profile (equation (6) with $`\alpha =1`$) for the virial mass $`M_v=10^{15}M_{\mathrm{}}`$ and with concentration $`c=6`$ (as suggested by formula (8) for the mass of this order). We see that for radial orbits the kurtosis profile is convex as opposed to the concave shapes in the case of isotropic and circular orbits. Since our data seem to prefer a concave shape we can expect isotropic or tangential orbits to fit the data best.

We begin by fitting the line-of-sight velocity dispersion profile shown in the upper panel of Figure 3. We consider different values of $`\alpha =0,1/2,1,3/2`$, and, for each of them, we determine the best-fitting anisotropy parameter $`\beta `$ and dark virial mass $`M_v`$ as a function of concentration $`c`$. The best-fitting values of $`M_v`$ obtained are of the order of $`10^{15}M_{\mathrm{}}`$, which corresponds to dark halo virial radius $`r_v=88^{}`$ or $`2.6`$ Mpc. The virialised region is supposed to lie within this radius, so in the following analysis we discard the two outer radial bins.

Figure 4 shows the results of fitting just the inner $`7`$ data points in the upper panel of Figure 3 for $`\alpha =0,1/2,1`$ and $`3/2`$. The best-fitting dark virial mass $`M_v`$ (upper panel) and the velocity anisotropy parameter for early-type galaxies $`\beta `$ (middle panel) are shown as a function of concentration. The lower panel of Figure 4 shows the goodness of fit $`\chi ^2`$ (we do not use $`\chi ^2/N`$, with $`N`$ the number of degrees of freedom, because it is not obvious how many parameters can be estimated with this procedure). As can be seen in Figure 4, the biggest $`M_v`$ obtained is $`1.5\times 10^{15}M_{\mathrm{}}`$ which gives $`r_v=100^{}`$ or $`2.9`$ Mpc. The data points are now all within this region and therefore there is no need for further adjustments in the number of data points analysed. The lowest panel of Figure 4 proves that neither $`c`$ nor $`\alpha `$ can be constrained from the analysis of velocity dispersion alone: $`\chi ^2`$ flattens for large $`c`$ and reaches similar values for all $`\alpha `$ for a wide range of $`c`$.

As discussed in Appendix A, the sampling distributions of $`\sigma _{\mathrm{los}}`$ and $`(\mathrm{log}\kappa _{\mathrm{los}})^{1/10}`$ are independent, hence we can use the same $`\chi ^2`$ minimisation scheme to find joint constraints following from fitting both quantities. Using the total of $`14`$ data points at $`R<80^{}`$, we jointly fit our four parameters, $`M_v`$, $`\beta `$, $`\alpha `$ and $`c`$. Contrary to the case when only velocity dispersion was studied, the minimisation procedure now converges. The minimum is found at $`M_v=1.2\times 10^{15}M_{\mathrm{}}`$ (corresponding to dark matter virial radius $`r_v=92^{}=2.7`$ Mpc), $`\beta =0.13`$, $`\alpha =0`$ and $`c=19`$ with $`\chi ^2/N=6.1/10`$. For a better visualisation of the constraints obtained for our 4 parameters, we plot in Figure 5 the cuts through the 4-dimensional confidence region in all six possible planes with probability contours corresponding to $`1\sigma `$ ($`68`$%), $`2\sigma `$ ($`95`$%) and $`3\sigma `$ ($`99.7`$%) i.e. $`\mathrm{\Delta }\chi ^2=\chi ^2\chi _{\mathrm{min}}^2=4.72,9.70,16.3`$, where $`\chi _{\mathrm{min}}^2=6.1`$.

Although the cuts do not tell everything about the confidence region, Figure 5 can be used to draw a number of interesting conclusions. The most striking is the behaviour of the confidence region in the $`c\alpha `$ plane shown in the upper right corner. Its shape shows that there is a strong degeneracy between the two parameters and indeed almost equally good fits can be obtained for the values of the inner slope other than $`\alpha =0`$. The best-fitting values of the remaining parameters (together with the corresponding $`\chi ^2`$ value) when different $`\alpha `$ are assumed are listed in Table 1. The results confirm that indeed $`\alpha =0`$ gives the best fit, but other inner slopes cannot be excluded. The steeper the inner slope (the higher the value of $`\alpha `$), however, the lower is the concentration required to provide good fits to the moments.

The velocity moments obtained with the sets of parameters listed in Table 1 are shown in Figure 6 together with the data; they overlap almost exactly. The dark matter profiles following from equation (6) with the parameters from Table 1 are plotted in the upper panel of Figure 7. In the lower panel we also show the logarithmic slopes of the profiles. In both panels our best-fitting profiles are compared to the “standard” NFW profile with concentration $`c=6`$ (as suggested by formula (8) for the mass of the order of $`10^{15}M_{\mathrm{}}`$). As can be seen in the upper panel, all our profiles overlap in a wide range of radial distances $`0.03r_v<r<r_v`$, and are somewhat steeper beyond $`0.05r_v`$ (see the lower panel) than the standard NFW model deduced from cosmological $`N`$-body simulations. Indeed, our best fitting NFW profile has concentration $`c=9.4`$, higher by 50% than the standard NFW profile for which $`c=6`$ according to equation (8).

The remaining parameters $`M_v`$ and $`\beta `$ are better constrained. The three panels on the left in Figure 5 show the probability contours for $`M_v`$ in the planes with the other 3 parameters. We find the best estimate for the dark mass to be $`M_v=(1.2\pm 0.4)\times 10^{15}M_{\mathrm{}}`$ ($`1\sigma `$ error-bars). The anisotropy parameter $`\beta `$ is very close to isotropic, the best fitting value $`\beta =0.13`$ gives $`\sigma _\theta =1.06\sigma _r`$, i.e. the best-fitting orbits of early-type galaxies are very weakly tangential, although fully consistent with isotropy, while radial orbits are excluded at the $`2\sigma `$ level.

The upper panel of Figure 8 compares the mass distributions in stars, gas and dark matter with the total mass distribution for our best fitting model with $`\alpha =0`$. At the dark matter virial radius $`r_v=2.7`$ Mpc, the total mass is $`1.4\times 10^{15}M_{\mathrm{}}`$. The mass in galaxies is only 2% of the total mass, the mass in gas is 13% of the total (the galaxies thus represent less than one-seventh of the baryonic mass), and the dark matter contributes the remaining 85% of the total mass. Therefore, the virial radius for the total mass is $`0.85^{1/3}=1.06`$ times that of the dark matter, i.e. $`97^{}`$ or $`2.9`$ Mpc.

The cumulative mass-to-light ratio, i.e. the ratio of the total mass distribution to the luminosity distribution in galaxies, is

$$M/L_B=\frac{M(r)}{L_\mathrm{G}(r)}.$$
(30)
where $`M(r)`$ is given by equation (29) with (2), (5) and (17), while $`L_\mathrm{G}(r)=M_\mathrm{G}(r)/\mathrm{{\rm Y}}`$. In the lower panel of Figure 8 we show $`M/L_B`$ for our best-fitting models for different $`\alpha `$, with parameters from Table 1. The models differ towards the centre and only for $`\alpha =1`$ does the distribution tend to a constant value in this limit, since we have used the NFW profile to fit the luminosity distribution of galaxies in Section 2. At radial distances larger than $`0.3`$ Mpc, the cumulative mass-to-light ratio decreases slowly to reach $`M/L_B351M_{\mathrm{}}/L_{\mathrm{}}`$ at the total mass virial radius $`2.9`$ Mpc.

## 6 Discussion

We studied the velocity moments of early-type galaxies in the Coma cluster and used them to constrain the distribution of dark matter and velocity anisotropy. Our analysis differs from previous analysis of optical data (e.g. The & White 1986; Merritt 1987; den Hartog & Katgert 1996; Carlberg et al. 1997; van der Marel et al. 2000; Biviano et al. 2003; Biviano & Girardi 2003), in that

1. we have, for a single cluster, a larger sample of galaxies, which, given their early morphological type, should be in dynamical equilibrium in the cluster potential;
2. we remove pairs from the computation of the velocity moments;
3. we include kurtosis in the analysis;
4. we model dark matter distribution using a generalised formula inspired by the results of cosmological $`N`$-body simulations;
5. we include hot gas in our analysis.

In comparison to studies based upon stacking of many clusters, our analysis of the Coma cluster benefits from not having to introduce errors in any stacking procedure, and from a cleaner removal of interlopers. On the other hand, the analyses of stacked clusters have the advantage of averaging out particular inhomogeneities of individual clusters such as Coma, expected in hierarchical scenarios of structure formation. Indeed, the Coma cluster is known to have irregular structure both in projected space (Fitchett & Webster 1987; Mellier et al. 1988; Briel et al. 1992) and velocity space (Colless & Dunn 1996; Biviano et al. 1996). In particular, the cluster has two central cD galaxies, NGC 4874 and NGC 4889, of which the first one is the central galaxy of the main cluster and the second probably belonged to a subcluster which has recently merged with the main cluster (e.g. Colless & Dunn 1996). There are other subgroups, such as the one associated with NGC 4839 at around $`40^{}`$ from NGC 4874, close enough to the cluster centre to contribute to our analysis.

The question whether the E-S0 sample is relaxed and how the existing substructure may affect our results can only be fully addressed by cosmological $`N`$-body simulations including galaxy formation, where all 3-dimensional information would be available. Although such an analysis has not yet been performed, the effect of the incomplete virialisation of structures of dark matter particles seen in cosmological simulations on the estimates of the mass of a single cluster through the Jeans equation has been addressed by Tormen, Bouchet & White (1997). They showed (see the bottom row of their Fig. 17) that even for significantly perturbed haloes the mass $`M(r)`$ at distances larger than 2% of the virial radius inferred by the proper Jeans analysis is within 30% (r.m.s.) of the true mass and departs from it by less than 20% (r.m.s.) for average or relaxed haloes. Since the dark matter distribution is known to possess more substructure than is observed in the galaxy distribution (cosmological simulations predict many more Milky Way satellites than are observed, see e.g. Moore et al. 1999) and since structures cease to grow sooner in flat universes with a cosmological constant, in comparison with analogous structures growing within the Einstein-de Sitter model (assumed in the simulations by Tormen et al.), they are more regular today (see Thomas et al. 1998) than shown in the first three columns of the bottom row of Figure 17 in Tormen et al. We therefore believe that the discrepancy between our derived mass and the true cluster mass due to substructure and departures from equilibrium is significantly smaller than the uncertainty due to sampling errors of the velocity moments.

Our results for the dark and total mass of the cluster are consistent with previous estimates. Using a combination of X-ray and optical data Hughes (1989) found for his preferred model a total mass within $`3.6h_{70}^1\mathrm{Mpc}`$ (where we used the notation $`H_0=70h_{70}\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$) to be $`(1.3\pm 0.2)\times 10^{15}h_{70}^1M_{\mathrm{}}`$, but a much wider range of masses if more general mass distributions were allowed. Briel et al. (1992) using ROSAT observations of Coma and assuming hydrostatic equilibrium of the gas derived a mass within the same radial distance to be $`(1.3\pm 0.4)\times 10^{15}h_{70}^1M_{\mathrm{}}`$. Also, by analysing the infall patterns around Coma, Geller, Diaferio & Kurtz (1999) find a mass of $`(1.7\pm 0.4)\times 10^{15}h_{70}^1M_{\mathrm{}}`$ within the same distance to the cluster centre as above. Extrapolating our results to this radial distance, we find an enclosed mass of $`(1.6\pm 0.5)\times 10^{15}h_{70}^1M_{\mathrm{}}`$, in good agreement with these three earlier estimates.

We find a strong degeneracy between the inner slope and the concentration of the dark matter profile, with many combinations of the two reproducing our velocity profiles almost equally well. In the range of inner slopes $`0\alpha 3/2`$ we find that the best-fitting models have the two parameters related almost linearly as $`c=199.6\alpha `$ with very little variation in the remaining fitting parameters, $`M_v`$ and $`\beta `$. The particular shape of the degeneracy between $`c`$ and $`\alpha `$ is due to the specific properties of the family of dark matter profiles used in this analysis, coupled with our lack of velocity data at radii smaller than the scale radius ($`r_s`$) of the dark matter. Using smaller radial bins would allow us to probe smaller radial distances, but at the expense of larger sampling errors in the velocity moments. As is clear from Figure 7, the best-fitting dark matter profiles obtained for different inner slopes are almost the same for a wide range of distances (larger than 3% of the virial radius) and this can only be achieved with the combination of the parameters as given in Table 1. Note that our best fit to the data for the flat dark matter density profile might mean that the data may be even better fitted with the unphysical model where the dark matter density profile rises with radius in the cluster centre, i.e. $`\alpha <0`$ in equation (6).

Our best-fitting NFW profile ($`\alpha =1`$) has a concentration $`c=9.4`$ (in agreement with the above formula), 50% higher than $`c=6`$ found in $`N`$-body simulations. It should be kept in mind, however, that the concentration parameters in the simulations are subject to substantial scatter and that the formula (8) (giving $`c=6`$) at masses of the order of $`10^{15}M_{\mathrm{}}`$ is actually an extrapolation of results obtained for lower mass haloes (see Bullock et al. 2001). However, the “standard” NFW model with $`c=6`$ is still within our $`1\sigma `$ confidence region in Figure 5.

In comparison, fits of the NFW profile to X-ray data of clusters, assumed isothermal, by Ettori & Fabian (1999), rescaled by Wu & Xue (2000), as well as by Sato et al. (2000), both yield $`c4`$ for the mass that we find for Coma within the virial radius (the latter after correction from $`\mathrm{\Delta }_\mathrm{c}=200`$ to $`\mathrm{\Delta }_\mathrm{c}=102`$), which is only within our $`2\sigma `$ confidence region. It may be that the assumption of isothermal hot gas causes a lower concentration parameter.

In an analysis similar in spirit to ours, Biviano et al. (2003) report $`c=4\pm 2`$ for their stacked ensemble of 59 ENACS clusters, fitted to an NFW model for the total mass density. Similarly, Biviano & Girardi (2003) show that NFW models with $`c=5.5`$ are consistent with a stacked ensemble of 43 2dFGRS clusters. Moreover, the ENACS and 2dFGRS clusters analysed in both studies are on average less massive (i.e. with lower velocity dispersions) than Coma, and given that cosmological simulations find a decreasing $`c`$ for increasing mass, the discrepancy with our result is even stronger.

The difference of these four studies with our result may be explained by the fact that the authors mentioned above fit the *total* mass density profile to the NFW form, while our fit was done for the distribution of the dark component only. If the gas is distributed similarly in the X-ray, ENACS and 2dFGRS clusters as in Coma, i.e. it has a flat inner core, and the gas on average contributes a substantial part of the total mass, fitting the total mass distribution may result in flatter profiles than the dark haloes really have. The discrepancy may be also caused by the exclusion of central cDs by Biviano et al. (2003), although we feel that removing one or two galaxies from our central bin of 39 galaxies will not affect our results. Alternatively, their mass model (obtained with the assumption of isotropic orbits) may not be coherent with the kurtosis profile of their stacked cluster.

Our best fit models have orbits that are very close to isotropic, as is expected, in the central regions, where the two-body relaxation time for the galaxy system is considerably smaller than the age of the Universe, and also because cosmological simulations indicate that dark matter particles typically have isotropic orbits in the centres of clusters (Thomas et al. 1998; Huss, Jain & Steinmetz 1999). In other words, we find no significant anisotropy bias for the galaxies relative to the expectations of isotropy for the dark matter. Note that if we force isotropic orbits for the galaxies, we then obtain similar constraints for the inner slope, concentration parameter and mass of the dark matter component within the virial radius as those shown in Figure 5.

One might wish to better reproduce the kurtosis profile. It is clear from Figure 6 that, while the velocity dispersion profile is well reproduced, there is room for improvement for the kurtosis profile. However, as discussed in Section 4, the kurtosis is mainly sensitive to the velocity anisotropy, which was modelled here by a single constant parameter. Besides, the curves shown in Figure 6 are the results of joint fitting of both moments and do not aim at reproducing the kurtosis alone. In spite of this, the inclusion of the kurtosis in our analysis allowed us to constrain the velocity anisotropy and other parameters of the model.

From our best estimate of the mass distribution, the baryons (galaxies and gas) contribute 15% of the total mass at the dark matter virial radius $`r_v=2.7`$ Mpc (see Figure 8). If we assume that the cluster content is representative of the Universe as a whole, we can use this baryon fraction to estimate the density parameter $`\mathrm{\Omega }_0`$ (see e.g. White et al. 1993). Taking the baryonic density parameter at its currently best value from nucleosynthesis $`\mathrm{\Omega }_b=0.04`$ (with $`h=0.7`$) we obtain $`\mathrm{\Omega }_0=0.26\pm 0.09`$ where the error comes from our 30% uncertainty in the dark virial mass value, 20% uncertainty in the gas mass (as estimated by White et al. 1993) and 10% error in the $`\mathrm{\Omega }_b`$ value (Burles, Nollett & Turner 2001).

Similarly, we may assume that clusters are good tracers, within their virial radius, of the ratio of mass to blue luminosity. Given that the closure mass-to-light ratio (critical density over luminosity density) in the blue band is roughly $`1100h_{70}`$ to 10% accuracy, consistent with the recent estimates of luminosity density from the 2dFGRS (Norberg et al. 2002) and SDSS (Blanton et al. 2001) galaxy surveys (after correction to the blue band), our mass-to-light ratio within the virial radius of $`351h_{70}`$ with 30% accuracy yields $`\mathrm{\Omega }_0=0.32\pm 0.1`$.

Combining these two estimates of the density parameter, we arrive at $`\mathrm{\Omega }_0=0.29\pm 0.1`$ in excellent agreement with other determinations, for example the recent value obtained from the WMAP CMB experiment (Spergel et al. 2003).

## Acknowledgements

We wish to thank Y. Hoffman and the anonymous referee for useful comments and suggestions. We acknowledge hospitality of the Institut d’Astrophysique de Paris where ELŁ benefited from the Jumelage exchange programme as well as the NATO Advanced Fellowship and where most of this work was done. Partial support was obtained from the Polish State Committee for Scientific Research within grant No. 2P03D02319. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We have also used the LEDA database (http://leda.univ-lyon1.fr) and SIMBAD operated at CDS, Strasbourg, France, as well as NASA’s Astrophysics Data System Service (ADS).

## Appendix A The sampling distributions

The validity of modelling presented in this paper rests on proper estimation of the velocity moments and their errors from observations. If their sampling distributions tend to normality, the statistics obtained from a sample (e.g. moments) can be characterised by an expectation value and a ‘standard’ error (Stuart & Ord 1994). For simplest statistics, like the second central moment, those can be calculated exactly from the population moments assuming the population properties. Since we are interested here also in quantities whose standard errors are not easily calculated analytically, we resort to Monte Carlo methods.

The most natural estimators of the variance and kurtosis from a sample of $`n`$ line-of-sight velocity measurements $`v_i`$ are

$$S^2=\frac{1}{n}\underset{i=1}{\overset{n}{}}(v_i\overline{v})^2$$
(31)
and

$$K=\frac{\frac{1}{n}_{i=1}^n(v_i\overline{v})^4}{(S^2)^2}$$
(32)
where

$$\overline{v}=\frac{1}{n}\underset{i=1}{\overset{n}{}}v_i$$
(33)
is the mean of galaxy velocities in the sample.

To investigate the distribution of these estimators for our binning of galaxies, i.e. when $`n=39`$, we ran Monte Carlo simulations by selecting $`𝒩=10^4`$ times $`n=39`$ numbers from a Gaussian distribution with zero mean and dispersion of unity. Since velocity distributions of gravitationally bound objects in general do not dramatically depart from a Gaussian, this is a sufficient approximation for constructing unbiased estimators of moments.

For each of the $`𝒩`$ samples we compute the statistic $`\theta _j^{}`$, namely $`S^2`$ or $`K`$ according to prescriptions given by equations (31) and (32). The Monte Carlo estimate of our statistic is then the mean of all values obtained

$$\theta ^{}=\frac{1}{𝒩}\underset{j=1}{\overset{𝒩}{}}\theta _j^{}$$
(34)
and its variance is

$$\mathrm{var}(\theta ^{})=\frac{1}{𝒩1}\underset{j=1}{\overset{𝒩}{}}(\theta _j^{}\theta ^{})^2.$$
(35)
We find that the best estimates obtained in this way are biased, especially for kurtosis (it is interesting to note, that the value of kurtosis is underestimated by a few percent even for much larger samples with $`n`$ of the order of few hundred). In addition, while the sampling distribution of velocity dispersion is Gaussian to a very good approximation, the one for kurtosis is strongly skewed. Using this knowledge we construct unbiased and Gaussian-distributed estimators of line-of-sight velocity dispersion $`s`$ and kurtosis-like variable $`k`$

$$s=\left(\frac{n}{n1}S^2\right)^{1/2}$$
(36)
$$k=\left[\mathrm{log}\left(\frac{3}{2.75}K\right)\right]^{1/10}.$$
(37)
The factor $`n1`$ in equation (36) is the well known correction for bias when estimating the sample variance, valid independently of the underlying distribution. In (37) the factor $`3/2.75`$ corrects for the bias in the kurtosis estimate, i.e. unbiased estimate of kurtosis is $`K^{}=3K/2.75`$, while the rather complicated function of $`K^{}`$ assures that the sampling distribution of $`k`$ is approximately Gaussian. Examples of sampling distributions of $`s`$ and $`k`$ from our Monte Carlo simulation are shown in Figure 9. We find that the standard errors in the case of $`s`$ are of the order of 11% (in agreement with an analytic result derivable with the assumption of Gaussian velocity distribution) while in the case of $`k`$ are approximately 2%.

The measured values of $`\sigma _{\mathrm{los}}`$ and $`\kappa _{\mathrm{los}}`$ calculated from our velocity data using equations (36)-(37) and (31)-(32) are shown in Figure 3. The $`1\sigma `$ error bars for the velocity dispersion are $`0.11s`$. The values of kurtosis are $`K^{}=3K/2.75`$ with approximate $`1\sigma `$ error bars propagated from the 2% error in $`k`$.

It is also important to check whether the sampling distributions of the two statistics are independent. In general, the covariance between, e.g. even moments derived from the same sample does not vanish, even for Gaussian distributions (Stuart & Ord 1994). However, the lowest order term is expected to decrease with the size of the sample $`n`$. To check whether $`n=39`$ in our case is indeed large enough to assure independence of $`s`$ and $`k`$, we construct the joint sampling distribution of the two statistics. The joint distribution in the form of $`𝒩`$ points with coordinates given by $`(s,k)`$ pairs calculated from each sample is presented in Figure 10. As indicated by the Figure the variables are very weakly correlated, we find the correlation coefficient of $`|\varrho |0.02`$.

To check the behaviour of $`s`$ and $`k`$ for velocity distributions departing from gaussianity, we repeated the Monte Carlo simulation, again sampling $`n=39`$ numbers from a Gaussian distribution, but modifying each sample by removing the 6 most inner points and adding 3 uniformly distributed in the range $`(1\sigma ,2\sigma )`$ and 3 in the range $`(2\sigma ,1\sigma )`$. Such distributions have unbiased kurtosis estimates of the order of $`K^{}=2.2`$, close to the lowest value obtained from our data (and most strongly departing from the Gaussian value of 3). We find that estimators $`s`$ and $`k`$ are again Gaussian-distributed to a very good approximation and very weakly correlated with $`|\varrho |0.07`$. The sampling distributions of $`s`$ and $`k`$ as well as the joint distribution look very similar to the purely Gaussian case. The correlation coefficient can increase significantly only in the presence of additional outliers (more numerous than predicted by the Gaussian distribution) in the range $`(2\sigma ,3\sigma )`$ and $`(3\sigma ,2\sigma )`$. We have checked that the number of galaxies with velocities in this range is 1 or 2 in each bin, in excellent agreement with the Gaussian prediction.

We can therefore assume that, to a good approximation, all our data points measuring velocity dispersion and kurtosis are independent, which justifies the use of standard $`\chi ^2`$ minimisation to fit the models to the data.


text:
# The 𝑁=2 Supersymmetric Heavenly Equation and Its Super-Hydrodynamical Reduction

## 1 Introduction

The so-called “Heavenly Equation” in $`(1+2)`$ dimensions was at first introduced to describe solutions of the complexified Einstein equations . It has received a lot of attention  for its remarkable integrability properties. Even if, to our knowledge, no attempt so far in the literature has been made to analyze the corresponding situation in the case of supergravity, the supersymmetric heavenly equation has been introduced and discussed in the context of superintegrable models. It appears  as the continuum limit of a system of equations known as “supersymmetric Toda lattice hierarchy”, introduced in . Such system of equations and their hidden $`N=2`$ supersymmetric structures have been vastly investigated . On the other hand, it has been recently pointed out that the dimensional reduction of the $`(1+2)`$-dimensional heavenly equation to a $`(1+1)`$-dimensional system is related to Benney-like integrable hierarchies of hydrodynamical equations . Another very recent reference discussing the hydrodynamical properties of the reduced heavenly equation is .

In this paper we introduce the $`N=2`$ supersymmetric heavenly equation, clarifying its integrability properties and its algebraic arising from the $`n\mathrm{}`$ limit of a class of $`N=2`$ supersymmetric Toda equations. Further, we investigate its $`(1+1)`$ dimensional reduction which provides the supersymmetric hydrodynamical equations extending the Benney-like hierarchy of reference . The integrability properties of such systems are induced by the integrability property of the $`N=2`$ superheavenly equation and its Lax representation discussed below. This Lax representation is not in the form of supersymmetric dispersionless Lax operator, as we discuss at length later.

The plan of the paper is as follows, in the next section we introduce the discretized $`N=2`$ Toda lattice hierarchy (as well as its continuum limit, the $`N=2`$ Superheavenly Equation), as a supersymmetric Toda system based on the $`sl(n|n+1)`$ superalgebra. It is worth to remark that these are the same superalgebras originally employed in the construction of . On the other hand, the superalgebras of the $`sl(n|n+1)`$ series admit a complex structure, allowing the construction of superToda systems based on $`N=2`$ superfields, according to the scheme of . For this reason, the next section results provide a generalization of those of . The supersymmetric Lax pairs providing integrability of the systems (for any given $`n`$ and in the limit $`n\mathrm{}`$) are explicitly constructed. They provide a dynamical formulation in the $`x_\pm `$ plane (without involving the extra-time direction $`\tau `$ of the superheavenly equation). In the following, we investigate the dimensional reduction of both the discretized and continuum systems from $`(1+2)`$ to $`(1+1)`$ dimensions. We obtain a supersymmetric system of equations with an interesting and subtle property. Unlike its purely bosonic subsector, the whole system involving fermions is not of hydrodynamical type. However, the same system, once expressed in terms of a supergeometry involving superfields and fermionic derivatives, satisfies a graded extension of the hydrodynamical property, which can be naturally named of “super-hydrodynamical type of equation”. With this expression we mean that there is a way to recast the given ($`1+1`$)-dimensional supersymmetric equations into a system of non-linear equations for superfields which involves only first-order derivations w.r.t. the supersymmetric fermionic derivatives. This super-hydrodynamical system furnishes the supersymmetrization of the system introduced in .

## 2 The $`N=2`$ Superheavenly equation.

The construction of a continuum limit (for $`n\mathrm{}`$) of a discretized superToda system requires a presentation of the system in terms of a specific Cartan matrix. The symmetric choice in  for the Cartan matrix of the superalgebra $`sl(n|n+1)`$ does not allow to do so. On the other hand , the Cartan matrix $`a_{ij}`$ of $`sl(n|n+1)`$ can be chosen to be antisymmetric with the only non-vanishing entries given by $`a_{ij}=\delta _{i,i+1}\delta _{i,i1}`$.

The Cartan generators $`H_i`$ and the fermionic simple roots $`F_{\pm i}`$ satisfy

$`[H_i,F_{\pm j}]`$ $`=`$ $`\pm a_{ij}F_{\pm j},`$
$`\{F_i,F_j\}`$ $`=`$ $`\delta _{ij}.`$ (2.1)
The continuum limit of the  construction could have been performed for any superalgebra admitting an $`n\mathrm{}`$ limit, such as $`sl(n|n)`$, etc. On the other hand, the superalgebras of the series $`sl(n|n+1)`$ are special because they admit a complex structure and therefore the possibility of defining an $`N=2`$ manifestly supersymmetric Toda system, following the prescription of . This is the content of the present section.

At first we introduce the $`N=2`$ fermionic derivatives $`D_\pm `$, $`\overline{D}_\pm `$, acting on the $`x_\pm `$ $`2D`$ spacetime ($`\theta _\pm `$ and $`\overline{\theta }_\pm `$ are Grassmann coordinates). The $`2D`$ spacetime can be either Euclidean ($`x_\pm =x\pm t`$) or Minkowskian ($`x_\pm =x\pm it`$).

We have

$`D_\pm `$ $`=`$ $`{\displaystyle \frac{}{\theta _\pm }}i\overline{\theta }_\pm _\pm ,`$
$`\overline{D}_\pm `$ $`=`$ $`{\displaystyle \frac{}{\overline{\theta }_\pm }}+i\theta _\pm _\pm .`$ (2.2)
They satisfy the anticommutator algebra

$`\{D_\pm ,\overline{D}_\pm \}`$ $`=`$ $`2i_\pm `$ (2.3)
and are vanishing otherwise.

Chiral ($`\mathrm{\Phi }`$) and antichiral ($`\overline{\mathrm{\Phi }}`$) $`N=2`$ superfields are respectively constrained to fulfill the conditions

$`\overline{D}_\pm \mathrm{\Phi }`$ $`=`$ $`0,`$
$`D_\pm \overline{\mathrm{\Phi }}`$ $`=`$ $`0.`$ (2.4)
Accordingly, a generic chiral superfield $`\mathrm{\Phi }`$ is expanded in its bosonic $`\phi `$, $`F`$ (the latter auxiliary) and fermionic component fields ($`\psi _+,\psi _{}`$) as

$`\mathrm{\Phi }(\widehat{x}_\pm ,\theta _\pm )`$ $`=`$ $`\phi +\theta _+\psi _++\theta _{}\psi _{}+\theta _+\theta _{}F,`$ (2.5)
with $`\phi `$, $`\psi _\pm `$ and $`F`$ evaluated in $`\widehat{x}_\pm =x_\pm +i\overline{\theta }_\pm \theta _\pm `$.

Similarly, the antichiral superfield $`\overline{\mathrm{\Phi }}`$ is expanded as

$`\overline{\mathrm{\Phi }}(\overline{x}_\pm ,\overline{\theta }_\pm )`$ $`=`$ $`\overline{\phi }+\overline{\theta }_+\overline{\psi }_++\overline{\theta }_{}\overline{\psi }_{}+\overline{\theta }_+\overline{\theta }_{}\overline{F},`$ (2.6)
with all component fields evaluated in $`\overline{x}_\pm =x_\pm i\overline{\theta }_\pm \theta _\pm `$.

Due to the complex structure of $`sl(n|n+1)`$, its Cartan and its simple (positive and negative) root sector can be split into its conjugated parts

$`\begin{array}{cc}\{H_{2k1}\},\hfill & \overline{}\{H_{2k}\},\hfill \\ _+\{F_{2k1}\},\hfill & _{}\{F_{(2k1)}\},\hfill \\ \overline{}_+\{F_{2k}\},\hfill & \overline{}_{}\{F_{2k}\},\hfill \end{array}`$ (2.10)
for $`k=1,2,\mathrm{},n`$ .

Following , we can introduce the $`sl(n|n+1)`$ $`N=2`$ superToda dynamics, defined for the Cartan-valued chiral ($`𝚽`$) and antichiral ($`\overline{𝚽}`$ ) $`N=2`$ superfields,

$`𝚽`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{n}{}}}\mathrm{\Phi }_kH_{2k1},`$
$`\overline{𝚽}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{n}{}}}\overline{\mathrm{\Phi }}_kH_{2k},`$ (2.11)
through the Lax operators $`_\pm `$ and $`\overline{}_\pm `$, given by

$`_+`$ $`=`$ $`D_+𝚽+e^{\overline{𝚽}}F_+e^{\overline{𝚽}},`$
$`_{}`$ $`=`$ $`F_{}`$ (2.12)
and

$`\overline{}_+`$ $`=`$ $`\overline{F}_+,`$
$`\overline{}_{}`$ $`=`$ $`\overline{D}_{}\overline{𝚽}+e^𝚽\overline{F}_{}e^𝚽,`$ (2.13)
where

$`\begin{array}{cc}F_+=_kF_{2k1},\hfill & F_{}=_kF_{(2k1)},\hfill \\ \overline{F}_+=_kF_{(2k)},\hfill & \overline{F}_{}=_kF_{2k},\hfill \end{array}`$ (2.16)
(as before, the sum is over the positive integers up to $`n`$).

Explicitly, we have

$`_+`$ $`=`$ $`{\displaystyle \underset{k}{}}(D_+\mathrm{\Phi }_kH_{2k1}+e^{\overline{\mathrm{\Phi }}_{k1}\overline{\mathrm{\Phi }}_k}F_{2k+1}),`$
$`_{}`$ $`=`$ $`{\displaystyle \underset{k}{}}F_{(2k1)},`$
$`\overline{}_+`$ $`=`$ $`{\displaystyle \underset{k}{}}F_{2k},`$
$`\overline{}_{}`$ $`=`$ $`{\displaystyle \underset{k}{}}(\overline{D}_{}\overline{\mathrm{\Phi }}_kH_{2k}+e^{\mathrm{\Phi }_k\mathrm{\Phi }_{k+1}}F_{2k}),`$
Please notice that here and in the following we have formally set $`\overline{\mathrm{\Phi }}_00`$.

The zero-curvature equations are given by

$`D_+_{}+D_{}_++\{_+,_{}\}`$ $`=`$ $`0,`$
$`\overline{D}_+\overline{}_{}+\overline{D}_{}\overline{}_++\{\overline{}_+,\overline{}_{}\}`$ $`=`$ $`0,`$ (2.18)
so that the following set of equations for the constrained (anti)chiral $`N=2`$ superfields is obtained

$`D_+D_{}\mathrm{\Phi }_k`$ $`=`$ $`e^{\overline{\mathrm{\Phi }}_{k1}\overline{\mathrm{\Phi }}_k},`$
$`\overline{D}_+\overline{D}_{}\overline{\mathrm{\Phi }}_k`$ $`=`$ $`e^{\mathrm{\Phi }_k\mathrm{\Phi }_{k+1}},`$ (2.19)
for the positive integers $`k=1,2,\mathrm{},n`$.

By setting,

$`B_k=\mathrm{\Phi }_k\mathrm{\Phi }_{k+1},`$ $`\overline{B}_k=\overline{\mathrm{\Phi }}_k\overline{\mathrm{\Phi }}_{k+1},`$ (2.20)
we get the two systems of equations

$`D_+D_{}B_k=e^{\overline{B}_k}e^{\overline{B}_{k1}},`$ $`\overline{D}_+\overline{D}_{}\overline{B}_k=e^{B_{k+1}}e^{B_k},`$ (2.21)
for $`k=1,2,\mathrm{},n`$.

By identifying $`k`$ as a discretized extra time-like variable $`\tau `$ we obtain, in the continuum limit for $`n\mathrm{}`$,

$`D_+D_{}B=_\tau e^{\overline{B}},`$ $`\overline{D}_+\overline{D}_{}\overline{B}=_\tau e^B,`$ (2.22)
which corresponds to the $`N=2`$ extension of the superheavenly equation.

Indeed, the presence in the previous equations of the first derivative in $`\tau `$ is an artifact of the $`N=2`$ superfield formalism. Once solved the equations at the level of the component fields and eliminated the auxiliary fields in terms of the equations of motion, we are left with a system of second-order equations.

We have, in terms of the component fields,

$`B`$ $`=`$ $`\left(1+i\overline{\theta }_+\theta _+_++i\overline{\theta }_{}\theta _{}_{}\overline{\theta }_+\theta _+\overline{\theta }_{}\theta _{}_+_{}\right)𝒞,`$
$`\overline{B}`$ $`=`$ $`\left(1i\overline{\theta }_+\theta _+_+i\overline{\theta }_{}\theta _{}_{}\overline{\theta }_+\theta _+\overline{\theta }_{}\theta _{}_+_{}\right)\overline{𝒞},`$ (2.23)
where

$`𝒞`$ $`=`$ $`\left(b+\theta _+\psi _++\theta _{}\psi _{}+\theta _+\theta _{}a\right),`$
$`\overline{𝒞}`$ $`=`$ $`\left(\overline{b}+\overline{\theta }_+\overline{\psi }_++\overline{\theta }_{}\overline{\psi }_{}+\overline{\theta }_+\overline{\theta }_{}\overline{a}\right),`$ (2.24)
with $`a`$, $`\overline{a}`$ bosonic auxiliary fields. All component fields are evaluated in $`x_\pm `$ only.

The equations of motion of the $`N=2`$ superheavenly equation are explicitly given in components by

$`a`$ $`=`$ $`(e^{\overline{b}})_\tau ,`$
$`2i_{}\psi _+`$ $`=`$ $`(\overline{\psi }_{}e^{\overline{b}})_\tau ,`$
$`2i_+\psi _{}`$ $`=`$ $`(\overline{\psi }_+e^{\overline{b}})_\tau ,`$
$`4_+_{}b`$ $`=`$ $`(\overline{a}e^{\overline{b}}\overline{\psi }_+\overline{\psi }_{}e^{\overline{b}})_\tau ,`$ (2.25)
and

$`\overline{a}`$ $`=`$ $`(e^b)_\tau ,`$
$`2i_{}\overline{\psi }_+`$ $`=`$ $`(\psi _{}e^b)_\tau ,`$
$`2i_+\overline{\psi }_{}`$ $`=`$ $`(\psi _+e^b)_\tau ,`$
$`4_+_{}\overline{b}`$ $`=`$ $`(ae^b\psi _+\psi _{}e^b)_\tau .`$ (2.26)
Eliminating the auxiliary fields we obtain the systems

$`2i_{}\psi _+`$ $`=`$ $`(\overline{\psi }_{}e^{\overline{b}})_\tau ,`$
$`2i_+\psi _{}`$ $`=`$ $`(\overline{\psi }_+e^{\overline{b}})_\tau ,`$
$`4_+_{}b`$ $`=`$ $`((e^b)_\tau e^{\overline{b}}+(\overline{\psi }_+\overline{\psi }_{}e^{\overline{b}})_\tau ,`$ (2.27)
and

$`2i_{}\overline{\psi }_+`$ $`=`$ $`(\psi _{}e^b)_\tau ,`$
$`2i_+\overline{\psi }_{}`$ $`=`$ $`(\psi _+e^b)_\tau ,`$
$`4_+_{}\overline{b}`$ $`=`$ $`((e^{\overline{b}})_\tau e^b+(\psi _+\psi _{}e^b)_\tau .`$ (2.28)
The bosonic component fields $`b`$, $`\overline{b}`$, as well as the fermionic ones $`\psi _\pm `$, $`\overline{\psi }_\pm `$, are all independent. The equations (2, 2) are a manifestly $`N=2`$ supersymmetric extension of the system introduced in .

## 3 Super-hydrodynamical reductions of the superheavenly equation.

The equations (2) and (2) correspond to a $`(1+2)`$-dimensional system, manifestly relativistic and $`N=2`$ supersymmetric w.r.t. the two-dimensional subspace spanned by the $`x_\pm `$ coordinates, while possessing an extra bosonic time-like dimension denoted as $`\tau `$. A very interesting example of integrable system in $`1+1`$-dimension which is currently intensely investigated, see e.g. , can be recovered by applying a dimensional-reduction to, let’s say, the bosonic sector of the $`(1+2)`$-dimensional heavenly equation. We can refer to such a system, perhaps a bit pedantically, as the $`(1+1)`$-dimensionally reduced heavenly equation. It can be obtained by setting equal to zero the fermionic variables $`\psi _\pm ,\overline{\psi }_\pm 0`$, while $`b`$, $`\overline{b}`$, can be consistently constrained as $`\overline{b}=b`$. The $`x_\pm `$ coordinates are identified, i.e. $`x_+=x_{}=x`$.

The resulting equation, by changing the normalizations (setting $`f=2b`$ and $`t=2\tau `$) can be conveniently written as

$$f_{tt}=(e^f)_{xx}$$
(3.29)
This equation has recently received a lot of attention in the literature. It is an example of a completely integrable, hydrodynamical-type equation. It admits a multihamiltonian structure and possesses an infinite number of conserved charges in involution. It has been recently shown that it can be recovered via a dispersionless Lax representation given by

$$L:=p^1\left(1+gp^2+hp^4\right)^{\frac{1}{2}},$$
(3.30)
with $`f,g`$ functions of $`x,t`$, while $`p`$ is the classical momentum. The equations of motion are read from

$$\frac{L}{t}=\frac{p}{2}\{L_0^2,L\}$$
(3.31)
where $`\{,\}`$ denotes the usual Poisson brackets and $`L_0^2`$, defined in , is explicitly given by $`L_0^2=p^2+g`$.

The equation (3.31) leads to

$`{\displaystyle \frac{h}{t}}`$ $`=`$ $`hg_x`$ (3.32)
$`{\displaystyle \frac{g}{t}}`$ $`=`$ $`h_x`$
The equation (3.29) is recovered after eliminating $`g`$ from the previous system and setting $`h=e^f`$. In the (3) form, the “$`(1+1)`$-dimensionally reduced heavenly equation” is shown to be a hydrodynamical type of equation.

An interesting question, whose solution as we will see is non-trivial, is whether the $`(1+1)`$-dimensional reduction (for $`x_+=x_{}=x`$) of the full $`N=2`$ supersymmetric heavenly equation is also of hydrodynamical type. The answer is subtle. The introduction of the fermionic fields $`\psi _\pm `$, $`\overline{\psi }_\pm `$, whose equations of motion are of first order in the extra-bosonic time, does not allow to represent the dimensionally reduced system from (2, 2) in hydrodynamical form, due to a mismatch with the second-order-equation satisfied by the bosonic fields $`b`$, $`\overline{b}`$. On the other hand, it is quite natural to expect that important properties are not spoiled by the supersymmetrization. This is indeed the case with the hydrodynamical reduction, when properly understood. The key issue here is the fact that the nice features of the supersymmetry are grasped when inserted in the proper context of a supergeometry, which must be expressed through the use of a superfield formalism. It is in this framework that a super-hydrodynamical reduction of the dimensionally reduced system from (2, 2) becomes possible. Indeed, for $`x_+=x_{}=x`$, we can write down the (2.22) system through the set of, first-order in the fermionic derivatives, superfield equations

$`𝒟_{}`$ $`=`$ $`N_\tau ,`$ (3.33)
$`𝒟_+N`$ $`=`$ $`e^{\overline{B}}`$
and

$`\overline{𝒟}_{}\overline{}`$ $`=`$ $`\overline{N}_\tau ,`$ (3.34)
$`\overline{𝒟}_+\overline{N}`$ $`=`$ $`e^B`$
with $`N`$ and $`\overline{N}`$ subsidiary fermionic superfields . Taking into account the following expansions

$`e^{\overline{B}}`$ $`=`$ $`(1i\overline{\theta }_+\theta _+_+i\overline{\theta }_{}\theta _{}_{}\overline{\theta }_+\theta _+\overline{\theta }_{}\theta _{}_+_{})e^{\overline{𝒞}}=`$
$`𝒟_+\theta _+\left(1i\overline{\theta _{}}\theta _{}_{}\right)e^{\overline{𝒞}}`$
$`e^B`$ $`=`$ $`(1+i\overline{\theta }_+\theta _+_++i\overline{\theta }_{}\theta _{}_{}\overline{\theta }_+\theta _+\overline{\theta }_{}\theta _{}_+_{})e^𝒞=`$ (3.35)
$`\overline{𝒟}_+\overline{\theta }_+\left(1+i\overline{\theta }_{}\theta _{}_{}\right)e^𝒞`$
we easily obtain the following solutions for $`N`$, $`\overline{N}`$

$`N`$ $`=`$ $`\left(\theta _+i\theta _+\overline{\theta }_{}\theta _{}_{}\right)e^{\overline{𝒞}}+𝒟_+\mathrm{\Omega }_1`$
$`\overline{N}`$ $`=`$ $`\left(\overline{\theta }_++i\overline{\theta }_+\overline{\theta }_{}\theta _{}_{}\right)e^𝒞+\overline{𝒟}_+\overline{\mathrm{\Omega }}_1`$ (3.36)
in terms of arbitrary $`\mathrm{\Omega }_1`$, $`\mathrm{\Omega }_2`$ bosonic superfunctions.

If we set, as we are free to choose,

$`\mathrm{\Omega }_1`$ $`=`$ $`\theta _+\left(\mathrm{\Psi }_{}2i\overline{\theta }_{}b_{}i\overline{\theta }_{}\theta _{}\mathrm{\Psi }_,\right)`$
$`\overline{\mathrm{\Omega }}_1`$ $`=`$ $`\overline{\theta }_+\left(\overline{\mathrm{\Psi }}_{}2i\theta _{}\overline{b}_{}+i\overline{\theta }_{}\theta _{}\overline{\mathrm{\Psi }}_,\right)`$ (3.37)
and substitute the values of $`N`$, $`\overline{N}`$ given in (3.36) back to (3.33) and (3.34), we obtain the equivalence of this system of equations w.r.t. the $`x_+=x_{}=x`$ dimensional reduction of the (2) and (2) equations. This proves the existence of a super-hydrodynamical reduction expressed in a superfield formalism. If we express in terms of the component fields this set of equations, as mentioned before, we are not obtaining a hydrodynamical type equation. This fact should not be regarded as a vicious feature of our system, but rather as a virtue of the supersymmetry. In several examples, this is one, the introduction of the super-formalism allows extending both the properties of the systems and the techniques used to investigate them, to cases for which the ordinary methods are of no help. In a somehow related area we can cite, e.g., the derivation of bosonic integrable hierarchies associated with the bosonic sector of super-Lie algebras . They are outside the classification based on and cannot be produced from ordinary Lie algebras alone.

Let us finally discuss the restriction from $`N=2N=1`$, i.e. down to the $`N=1`$ supersymmetry. It is recovered from setting

$`\overline{\theta }_+=\theta _+,`$ $`\overline{\theta }_{}=\theta _{},`$
$`\overline{𝒟}_+=𝒟_+,`$ $`\overline{𝒟}_{}=𝒟_{}.`$ (3.38)
The equations (2) and (2) now read

$`𝒟_+𝒟_{}𝒞=_\tau e^{\overline{𝒞}},`$ $`𝒟_+𝒟_{}\overline{𝒞}=_\tau e^𝒞,`$ (3.39)
with

$$D_\pm =\frac{}{\theta _\pm }+i\theta _\pm _\pm $$
(3.40)
The equations in (3.39) have a similar structure in components as the equations (2) and (2). We have indeed

$`_{}\psi _+`$ $`=`$ $`i(\overline{\psi }_{}e^{\overline{b}})_\tau ,`$
$`_+\psi _{}`$ $`=`$ $`i(\overline{\psi }_+e^{\overline{b}})_\tau ,`$
$`_+_{}b`$ $`=`$ $`((e^b)_\tau e^{\overline{b}}+(\overline{\psi }_+\overline{\psi }_{}e^{\overline{b}})_\tau ,`$ (3.41)
and

$`_{}\overline{\psi }_+`$ $`=`$ $`i(\psi _{}e^b)_\tau ,`$
$`_+\overline{\psi }_{}`$ $`=`$ $`i(\psi _+e^b)_\tau ,`$
$`_+_{}\overline{b}`$ $`=`$ $`((e^{\overline{b}})_\tau e^b+(\psi _+\psi _{}e^b)_\tau .`$ (3.42)
If we further dimensionally reduce to $`(1+1)`$-dimension our $`N=1`$ system, by constraining the variables $`x_\pm `$ ($`x_+=x_{}=it`$) and setting $`\overline{𝒞}=𝒞`$, we obtain the system

$`_t\psi _+`$ $`=`$ $`_\tau \left(e^b\psi _{}\right)`$
$`_t\psi _{}`$ $`=`$ $`_\tau \left(e^b\psi _+\right)`$
$`_t^2b`$ $`=`$ $`_\tau \left({\displaystyle \frac{1}{2}}(e^{2b})_\tau +e^b\psi _+\psi _{}\right).`$ (3.43)
As its $`N=2`$ counterpart, this system is not of hydrodynamical type. However, it can be easily shown, with a straightforward modification of the procedure previously discussed for the $`N=2`$ case, to be super-hydrodynamical when expressed in terms of superfields and superderivatives.

## 4 Conclusions.

In this paper we have at first introduced manifest $`N=2`$ supersymmetric discrete Toda equations in association with the complex-structure superalgebras of the $`sl(n,n+1)`$ series. We proved that in the $`n\mathrm{}`$ limit the discrete variable can be regarded as a continuum extra bosonic time $`\tau `$. The corresponding system is a manifestly $`N=2`$ supersymmetric extension of the standard heavenly equation in $`(2+1)`$ dimensions, the supersymmetry being associated with the plan spanned by the bosonic $`x_\pm `$ coordinates and their fermionic counterparts, which are naturally expressed in an $`N=2`$ superspace formalism. A supersymmetric Lax-pair, proving the integrability of the system, was constructed. It expresses the dynamics w.r.t. the $`\pm `$ superspace coordinates ($`\tau `$ in this respect can be considered either as an auxiliary parameter or as the discrete label of the original formulation).

We further investigated the dimensional reduction of the above systems from the $`(2+1)`$-dimensional to the $`(1+1)`$-dimensional case. In the purely bosonic case the reduced systems of equations are of a special type, they are hydrodynamical systems of non-linear equations. The supersymmetric case is much subtler. This system is not of hydrodynamical type since the presence of the fermions spoils this property. However, when properly expressed in terms of a supergeometry (essentially, superfields and fermionic derivatives) is of super-hydrodynamical type (the precise meaning of this term has been explained in the Introduction). More than a vice, this can be considered as a virtue of the supersymmetry. It allows extending the notion of hydrodynamical equations beyond the realm of the systems ordinary allowed.

Finally, it is worth mentioning a formal, however challenging problem, concerning the supersymmetrization. While the integrability of the systems under consideration is automatically guaranteed by the given supersymmetric and relativistic Lax-pairs in the $`\pm `$ plane mentioned above, one can wonder whether a dispersionless Lax operator, directly expressing the $`\tau `$ dynamics (i.e., in terms of the bosonic extra-time) could be found. The answer is indeed positive for the purely bosonic sector. In the supersymmetric case, however, no closed dispersionless Lax operator is known at present. It is repeated here an analogous situation already encountered for the polytropic gas systems . The problem of constructing supersymmetric dispersionless Lax operators for these related systems is still open.

Acknowledgments.

Z.P. is grateful for the hospitality at CBPF, where this work was initiated, while F.T. is grateful for the hospitality at the Institute of Theoretical Physics of the University of Wrocław, where the paper has been finished.