format: text

text:
# Untitled Document

An Ill Posed Cauchy Problem

for a Hyperbolic System in Two Space Dimensions

Alberto Bressan

S.I.S.S.A. - Via Beirut 4, Trieste 34014, Italy

e-mail: bressan@sissa.it

1 - Introduction The theory of weak solutions for nonlinear conservation laws is now well developed in the case of scalar equations  and for one-dimensional hyperbolic systems . For systems in several space dimensions, however, even the global existence of solutions to the Cauchy problem remains a challenging open question. In this note we construct a conterexample showing that, even for a simple class of hyperbolic systems, in two space dimensions the Cauchy problem can be ill posed. The systems on $`IR^m`$ that we consider take the special form

$$\frac{}{t}u_i+\underset{\alpha =1}{\overset{m}{}}\frac{}{x_\alpha }\left(f_\alpha \left(|u|\right)u_i\right)=0i=1,\mathrm{},n.$$
$`(1.1)`$
For physical motivations, see  or . Given a sufficiently regular initial condition

$$u(0,x)=\overline{u}(x),$$
$`(1.2)`$
the solution of the Cauchy problem can be constructed as follows. 1. Set $`\rho |u|`$ and solve the Cauchy problem for a scalar conservation law

$$\rho _t+\underset{\alpha =1}{\overset{m}{}}\left(f_\alpha (\rho )\rho \right)_{x_\alpha }=0,\rho (0,x)=\left|\overline{u}(x)\right|.$$
$`(1.3)`$
This will provide the absolute value of the solution. 2. To find the angular component $`\theta u/|u|`$, consider the O.D.E.

$$\dot{x}=f\left(\rho (t,x)\right)$$
$`(1.4)`$
with $`f=(f_1,\mathrm{},f_m)`$. Assuming that (1.4) has a unique solution for almost every initial data $`x(0)=y`$, denote by $`tx(t)\mathrm{\Phi }_ty`$ the corresponding trajectory. 3. Given $`t>0`$ and $`xIR^m`$, call $`y\mathrm{\Phi }^tx`$ be the (unique) point such that $`\mathrm{\Phi }_ty=x`$. Assuming that the measurable maps $`\mathrm{\Phi }^t`$ are well defined (up to sets of zero measure), the angular component of the solution of (1.1)-(1.2) is then obtained as

$$\theta (t,x)=\overline{\theta }\left(\mathrm{\Phi }^tx\right),$$
$`(1.5)`$
where $`\overline{\theta }\overline{u}/|\overline{u}|`$.

If the solution $`\rho =\rho (t,x)`$ of the scalar conservation law (1.3) remains smooth (or piecewise smooth), it is easy to check that the above procedure actually yields a weak solution $`u`$ to the Cauchy problem (1.1)-(1.2). In general however, even if the initial data $`\rho (0,)`$ is smooth, regularity may be lost in finite time. It is thus interesting to understand what kind of assumptions can provide the existence of entropy weak solutions. A natural set of conditions is the following:

(A1) The flux function $`F(\rho )f(\rho )\rho `$ is Lipschitz continuous.

(A2) The initial data is bounded, measurable, and bounded away from zero, i.e.

$$\overline{u}𝐋^{\mathrm{}},0<a\left|\overline{u}(x)\right|b\text{for a.e. }xIR^m.$$

In the case of one space dimension, the above assumptions yield the global existence of a unique entropy weak solution to the Cauchy problem. Indeed, by the fundamental theorem of Kruzhkov , the equation (1.3) admits a unique entropy weak solution $`\rho =\rho (t,x)`$, satisfying

$$a\rho (t,x)b.$$
$`(1.6)`$
Let $`x_1(t)<x_2(t)`$ be any two solutions of (1.4). For every $`t>0`$ the conservation equation implies

$$_{x_1(t)}^{x_2(t)}\rho (t,x)𝑑x=_{x_1(0)}^{x_2(0)}\rho (0,x)𝑑x.$$

By (1.6) this yields

$$\frac{a}{b}\left(x_2(0)x_1(0)\right)x_2(t)x_1(t)\frac{b}{a}\left(x_2(0)x_1(0)\right).$$

Hence the maps

$$(t,y)\mathrm{\Phi }_ty,(t,x)\mathrm{\Phi }^tx$$

are both well defined and Lipschitz continuous. It is now clear that (1.5) yields the desired solution. The example constructed below will show that, under exactly the same assumptions (A1)-(A2) which provide the existence and uniqueness of solutions in the one-dimensional case, in two space dimensions the Cauchy problem (1.1)-(1.2) is not well posed. The reason is that, although the scalar conservation law (1.3) always admits a unique entropy solution, the corresponding O.D.E. in (1.4) has discontinuous right hand side. For a given time $`t>0`$, the flux function $`y\mathrm{\Phi }_ty`$ may not admit a well defined inverse, because of oscillation phenomena.

2 - A counterexample Call $`𝐞_1=(1,0)`$, $`𝐞_2=(0,1)`$ the canonical basis on $`IR^2`$. Consider the Lipschitz (piecewise affine) flux function

$$F(\rho )\{\begin{array}{cc}0\hfill & \text{if }\rho 1\text{,}\hfill \\ (1\rho )𝐞_1\hfill & \text{if }1\rho 2\text{,}\hfill \\ (\rho 3)𝐞_1\hfill & \text{if }2\rho 3\text{,}\hfill \\ (\rho 3)𝐞_2\hfill & \text{if }3\rho \text{.}\hfill \end{array}$$
$`(2.1)`$
Observe that

$$F^{}(\rho )=\{\begin{array}{cc}0\hfill & \text{if }\rho 1\text{,}\hfill \\ 𝐞_1\hfill & \text{if }1\rho 2\text{,}\hfill \\ 𝐞_1\hfill & \text{if }2\rho 3\text{,}\hfill \\ 𝐞_2\hfill & \text{if }3\rho \text{.}\hfill \end{array}$$
$`(2.2)`$
We can write $`F(\rho )=f(\rho )\rho `$, with

$$f(\rho )\{\begin{array}{cc}0\hfill & \text{if }\rho 1\text{,}\hfill \\ \rho ^1(1\rho )𝐞_1\hfill & \text{if }1\rho 2\text{,}\hfill \\ \rho ^1(\rho 3)𝐞_1\hfill & \text{if }2\rho 3\text{,}\hfill \\ \rho ^1(\rho 3)𝐞_2\hfill & \text{if }3\rho \text{.}\hfill \end{array}$$
$`(2.3)`$
In particular,

$$f(2)=\frac{1}{2}𝐞_1,f(3)=0,f(4)=\frac{1}{4}𝐞_2.$$
$`(2.4)`$
We now look at the O.D.E.

$$\dot{x}=f\left(\rho (t,x)\right)$$
$`(2.5)`$
in connection with special types of solutions $`\rho =\rho (t,x)`$ of the conservation law

$$\rho _t+_xF(\rho )=0.$$
$`(2.6)`$
For a given set $`QIR^2`$, we consider three entropy weak solutions of (2.6), namely

$$\rho ^{\mathrm{}}(t,x)3,$$

$$\rho ^{\mathrm{}}(t,x)\{\begin{array}{cc}4\hfill & \text{if }xt𝐞_2Q\text{,}\hfill \\ 3\hfill & \text{otherwise,}\hfill \end{array}$$
$`(2.7)`$
$$\rho ^{\mathrm{}}(t,x)\{\begin{array}{cc}2\hfill & \text{if }xt𝐞_1Q\text{,}\hfill \\ 3\hfill & \text{otherwise.}\hfill \end{array}$$
$`(2.8)`$
figure 1

We begin by studying what happens to the corresponding trajectories of (2.5) in the special case where $`Q`$ is a rectangle:

$$Q[0,a]\times [0,b]IR^2.$$

1. In the case $`\rho =\rho ^{\mathrm{}}`$, we trivially have

$$\dot{x}=0\mathrm{\Phi }_t^{\mathrm{}}y=y.$$

2. When $`\rho =\rho ^{\mathrm{}}`$, we have (fig. 1)

$$\dot{x}=\{\begin{array}{cc}𝐞_2/4\hfill & \text{if }x_1[0,a]\text{}x_2t[0,b]\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

Points on the vertical strip

$$\mathrm{\Gamma }_a\left\{(x_1,x_2);x_1[0,a]\right\}$$

are eventually shifted to the right by an amount $`b𝐞_1/3`$, while all other points do not move. More precisely:

$$\mathrm{\Phi }_t^{\mathrm{}}yy=\{\begin{array}{cc}b𝐞_1/3\hfill & \text{if }y_1[0,a],y_2b,t(y_2b)+4b/3,\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

This is illustrated in fig. 1. One has $`\rho =4`$ on $`Q`$ and $`\rho =3`$ outside. The rectangle $`Q`$ moves upward with unit speed. At any given time $`t`$, points which lie inside the rectangle move with speed $`\dot{x}=𝐞_2/4`$. The points which initially lie on the segment $`\gamma `$ are eventually displaced onto the set $`\mathrm{\Phi }^{\mathrm{}}\gamma `$. 3. When $`\rho =\rho ^{\mathrm{}}`$, we have

$$\dot{x}=\{\begin{array}{cc}𝐞_1/2\hfill & \text{if }x_1t[0,a]\text{}x_2[0,b]\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

Points on the horizontal strip

$$\mathrm{\Gamma }_b\left\{(x_1,x_2);x_2[0,b]\right\}$$

are eventually shifted to the left by an amount $`a𝐞_1/3`$, while all other points do not move. More precisely:

$$\mathrm{\Phi }_t^{\mathrm{}}yy=\{\begin{array}{cc}a𝐞_1/3\hfill & \text{if }y_2[0,b],y_1a,t(y_1a)+2a/3,\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

Starting with these elementary solutions, we can construct more complex ones. For example, in (2.7) can take

$$Q=Q^{\mathrm{}}\underset{k\mathrm{even}}{}[k,k+1]\times [k,k+3]$$

In the corresponding equation (2.5), this generates the shift

$$\mathrm{\Phi }^{\mathrm{}}yy=\{\begin{array}{cc}𝐞_2\hfill & \text{if }[[y_1]]\text{ is even,}\hfill \\ 0\hfill & \text{if }[[y_1]]\text{ is odd.}\hfill \end{array}$$

On the other hand, in (2.8) we can take

$$Q=Q^{\mathrm{}}\underset{k\mathrm{even}}{}[k,k+3]\times [k,k+1]$$

In the corresponding equation (2.5), this generates the shift

$$\mathrm{\Phi }^{\mathrm{}}yy=\{\begin{array}{cc}𝐞_1\hfill & \text{if }[[y_2]]\text{ is even,}\hfill \\ 0\hfill & \text{if }[[y_2]]\text{ is odd,}\hfill \end{array}$$

where $`[[s]]`$ denotes the largest integer $`s`$.

figure 2

In the following, we write the coordinates of a point $`x=(x_1,x_2)IR^2`$ in binary digits, say

$$x_1=\underset{k}{}\alpha _k2^k,x_2=\underset{k}{}\beta _k2^k.$$
$`(2.9)`$
with $`\alpha _k,\beta _k\{0,1\}`$. For every $`k`$, the previous analysis shows that one can construct sets $`Q_k^{\mathrm{}}`$, $`Q_k^{\mathrm{}}`$, $`\stackrel{~}{Q}_k^{\mathrm{}}`$ contained inside disjoint strips

$$\begin{array}{cc}\hfill Q_k^{\mathrm{}}& \left\{(x_1,x_2);x_2+x_12^k[28,\mathrm{\hspace{0.17em}32}]\right\},\hfill \\ & \\ \hfill Q_k^{\mathrm{}}& \left\{(x_1,x_2);x_2+x_12^k[22,\mathrm{\hspace{0.17em}27}]\right\},\hfill \\ & \\ \hfill \stackrel{~}{Q}_k^{\mathrm{}}& \left\{(x_1,x_2);x_2+x_12^k[16,\mathrm{\hspace{0.17em}21}]\right\},\hfill \end{array}$$
$`(2.10)`$
and such that the following holds. (i) Given the solution

$$\rho _k^{\mathrm{}}(t,x)\{\begin{array}{cc}4\hfill & \text{if }xt𝐞_2Q_k^{\mathrm{}}\text{,}\hfill \\ 3\hfill & \text{otherwise,}\hfill \end{array}$$

the corresponding equation (2.5) generates the shift

$$\mathrm{\Phi }^{\mathrm{}}yy=\{\begin{array}{cc}2^k𝐞_2\hfill & \text{if }\alpha _k=0\text{,}\hfill \\ 0\hfill & \text{if }\alpha _k=1\text{.}\hfill \end{array}$$

(ii) Given the solution

$$\rho _k^{\mathrm{}}(t,x)\{\begin{array}{cc}2\hfill & \text{if }x+t𝐞_1Q_k^{\mathrm{}}\text{,}\hfill \\ 3\hfill & \text{otherwise,}\hfill \end{array}$$

the corresponding equation (2.5) generates the shift

$$\mathrm{\Phi }^{\mathrm{}}yy=\{\begin{array}{cc}2^k𝐞_1\hfill & \text{if }\beta _k\beta _{k+1}\text{,}\hfill \\ 0\hfill & \text{if }\beta _k=\beta _{k+1}\text{.}\hfill \end{array}$$

(iii) Given the solution

$$\stackrel{~}{\rho }_k^{\mathrm{}}(t,x)\{\begin{array}{cc}4\hfill & \text{if }xt𝐞_2\stackrel{~}{Q}_k^{\mathrm{}}\text{,}\hfill \\ 3\hfill & \text{otherwise,}\hfill \end{array}$$

the corresponding equation (2.5) generates the shift

$$\stackrel{~}{\mathrm{\Phi }}^{\mathrm{}}yy=\{\begin{array}{cc}2^{k1}𝐞_2\hfill & \text{if }\alpha _k=0\text{,}\hfill \\ 0\hfill & \text{if }\alpha _k=1\text{.}\hfill \end{array}$$

As shown in fig. 2, the set $`Q_k^{\mathrm{}}`$ consists of an array of (gray) rectangles of size $`2^k\times 32^k`$, the set $`Q_k^{\mathrm{}}`$ consists of an array of (black) rectangles of size $`2^{1k}\times 2^k`$, while the set $`\stackrel{~}{Q}_k^{\mathrm{}}`$ consists of an array of (gray) rectangles of size $`2^k\times 32^{1k}`$. Moving the gray rectangles upward with unit speed and moving the black ones to the right with unit speed, the three arrays will never overlap. We can thus consider the composite solution

$$\rho _k(t,x)\{\begin{array}{cc}4\hfill & \text{if }xt𝐞_2Q_k^{\mathrm{}}\stackrel{~}{Q}_k^{\mathrm{}}\text{,}\hfill \\ 2\hfill & \text{if }xt𝐞_1Q_k^{\mathrm{}}\text{,}\hfill \\ 3\hfill & \text{otherwise.}\hfill \end{array}$$
$`(2.11)`$
The differential equation

$$\dot{x}=f\left(\rho _k(t,x)\right)$$

will eventually shift points on the plane, according to the composition

$$\mathrm{\Psi }_k\stackrel{~}{\mathrm{\Phi }}_k^{\mathrm{}}\mathrm{\Phi }_k^{\mathrm{}}\mathrm{\Phi }_k^{\mathrm{}}.$$
$`(2.12)`$
figure 3

This is illustrated in fig. 3. Each square here has sides of length $`2^{k1}`$. In all four pictures, the asterisque marks the samepoint in the plane. Performing first the upward shift $`\mathrm{\Phi }_k^{\mathrm{}}`$, then the leftward shift $`\mathrm{\Phi }_k^{\mathrm{}}`$ and finally the upward shift $`\stackrel{~}{\mathrm{\Phi }}_k^{\mathrm{}}`$, one obtains a measure-preserving transformation $`\mathrm{\Psi }_k`$ whose main property is the following. In connection with the dyadic decomposition (2.9), define

$$X_j\left\{(x_1,x_2);\beta _j=0\right\},X_j^{}\left\{(x_1,x_2);\beta _j=1\right\}.$$

Then, up to sets of measure zero, one has

$$\mathrm{\Psi }_k(X_k)=X_{k+1},\mathrm{\Psi }_k(X_k^{})=X_{k+1}^{}.$$
$`(2.13)`$
Intuitively, one can think at the moving sets $`Q_k^{\mathrm{}}`$, $`Q_k^{\mathrm{}}`$, $`\stackrel{~}{Q}_k^{\mathrm{}}`$ as a comb, moving in the north-east direction. The effect of its passage is to displace points on the set $`X_k`$ onto the set $`X_{k+1}`$. Putting an arbitrarily large number of these combs one next to the other (fig. 4), we now construct a sequence of solutions $`u_n`$ where the initial data converge strongly in $`𝐋_{\mathrm{loc}}^1`$ but at a later time $`\tau `$ the corresponding solutions have no strong limit. Hence the limiting Cauchy problem cannot be well posed. Define a solution of (1.3) by setting

$$\rho _n(t,x)\{\begin{array}{cc}4\hfill & \text{if }xt𝐞_2_{0kn}\left(Q_k^{\mathrm{}}\stackrel{~}{Q}_k^{\mathrm{}}\right)\text{,}\hfill \\ & \\ 2\hfill & \text{if }xt𝐞_1_{0kn}Q_k^{\mathrm{}}\text{,}\hfill \\ & \\ 3\hfill & \text{otherwise.}\hfill \end{array}$$
$`(2.14)`$
Consider the Cauchy problem (1.1)-(1.2) with $`f`$ given at (2.3). We choose an initial value $`\overline{u}_n`$ of the form

$$\overline{u}_n(x)=\rho _n(0,x)(\mathrm{cos}\overline{\theta }(x),\mathrm{sin}\overline{\theta }(x))$$
$`(2.15)`$
with

$$\overline{\theta }(x)\{\begin{array}{cc}\beta \hfill & \text{if }[[x_2]]\text{ is even,}\hfill \\ \beta \hfill & \text{if }[[x_2]]\text{ is odd,}\hfill \end{array}$$

for a fixed angle $`\beta ]0,\pi /2[`$. We claim that the sequence of weak solutions $`u_n`$ has no limit in $`𝐋_{\mathrm{loc}}^1`$ as $`n\mathrm{}`$. Indeed, for $`32<x_1+x_2<t`$ one has

$$u_n(t,x)=3(\mathrm{cos}\theta _n(t,x),\mathrm{sin}\theta _n(t,x)),$$

with $`x=(x_1,x_2)`$ and

$$\theta _n(t,x)=\{\begin{array}{cc}\beta \hfill & \text{if }[[2^nx_2]]\text{ is even,}\hfill \\ \beta \hfill & \text{if }[[2^nx_2]]\text{ is odd.}\hfill \end{array}$$

As a result, when $`n\mathrm{}`$, the $`u_n`$ do not converge in $`𝐋_{\mathrm{loc}}^1`$ and the Cauchy problem is not well posed.

figure 4

On the other hand, still in the range where $`32<x_1+x_2<t`$, one can compute the weak limit of the sequence $`\{u_n\}`$, given by

$$u(t,x)=(3\mathrm{cos}\beta ,0).$$

To check whether this weak limit provides a distributions solution to the system (1.1), (2.3), it suffices to examine what happens in a neighborhood of the line $`x_1+x_2=t`$. Computing the time derivative

$$\frac{d}{dt}_{\mathrm{\Omega }(t)}u(t,x)𝑑x$$

for a set of the form

$$\mathrm{\Omega }(t)\{(x_1,x_2),x_2[a,b],x_1[tx_22^k,tx_2+2^k]\},$$

one checks that the conservation equations are satisfied provided that $`F(3\mathrm{cos}\beta )=F(3)=0`$. In other words: \- If $`\mathrm{cos}\beta 1/3`$, then the weak limit $`u`$ provides a weak solution to the system of conservation laws. This yields a curious example where entropy is dissipated by a linearly degenerate field. \- If $`\mathrm{cos}\beta >1/3`$, then the weak limit $`u`$ is not weak solution. In this case, it is not even clear whether any solution exists at all.

3 - Concluding remarks 1. All systems of the form (1.1) are hyperbolic, but in a certain way pathological. In one space dimension, the system

$$u_t+\left(f\left(|u|\right)u\right)_x=0$$
$`(3.1)`$
takes the quasilinear form

$$u_t+A(u)u_x=0,A(u)=f\left(|u|\right)I+f^{}\left(|u|\right)uu.$$

At a given point $`u0`$, the Jacobian matrix $`A(u)`$ has the 1-dimensional eigenspace $`EIRu`$ spanned by the vector $`u`$, corresponding to the eigenvalue $`\lambda =f\left(|u|\right)+f^{}\left(|u|\right)|u|`$. Moreover, one finds a second eigenvalue $`\lambda ^{}f\left(|u|\right)`$ of multiplicity $`n1`$, corresponding to the orthogonal eigenspace $`E^{}`$ (fig. 5).

figure 5

Even in one space dimension the Riemann problem can have multiple solutions, all obtained as vanishing viscosity limits. For example, consider

$$u_t+\left(|u^2|u\right)_x=0uIR^2,$$
$`(3.2)`$
$$u(0,x)=\{\begin{array}{cc}𝐞_1\hfill & \text{if }x<0\text{,}\hfill \\ 𝐞_1\hfill & \text{if }x>0\text{.}\hfill \end{array}$$

An entropic solution is provided by

$$u(t,x)=\rho (t,x)𝐞_1,$$

where $`\rho `$ is the entropy weak solution of the corresponding scalar problem

$$\rho _t+(\rho ^3)_x=0\rho (0,x)=\{\begin{array}{cc}1\hfill & \text{if }x<0\text{,}\hfill \\ 1\hfill & \text{if }x>0\text{.}\hfill \end{array}$$

This is clearly a limit of vanishing viscosity approximations .

In addition, there is the second admissible solution

$$u(t,x)=\{\begin{array}{cc}𝐞_1\hfill & \text{if }x<t\text{,}\hfill \\ 𝐞_1\hfill & \text{if }x>t\text{.}\hfill \end{array}$$
$`(3.3)`$
Indeed, consider any smooth function $`\theta :IR[0,\pi ]`$ such that

$$\theta (s)=\{\begin{array}{cc}0\hfill & \text{if }x0\text{,}\hfill \\ \pi \hfill & \text{if }x1\text{.}\hfill \end{array}$$

Then the Cauchy problem for (3.2) with initial data

$$u_n(0,x)=(\mathrm{cos}\theta (nx),\mathrm{sin}\theta (nx))$$

has a unique smooth solution, namely $`u_n(t,x)=u(0,xt)`$. This is also a limit of vanishing viscosity approximations. Letting $`n\mathrm{}`$, these solutions converge to (3.3).

We remark, however, that our counterexample is not related to the presence of a singularity at the origin. Indeed, the initial data at (2.14)-(2.15) are contained in a bounded domain whose convex closure does not contain the origin.

2. The flux function is not smooth but only Lipschitz continuous. The oscillations in the solution of the scalar conservation law (1.3) are not damped, because of the linear degeneracy of the flux. This indeed leaves open the possibility that, with a smooth flux, such a pathological behavior will not be observed.

We recall that, in the one dimensional case, it is common to replace a general flux $`F`$ by a polygonal approximation $`F_\nu `$. This is indeed the basis for the method of front-tracking approximations. In one space dimension, the Cauchy problems related to all polygonal fluxes are well posed . Taking a limit, one thus obtains a proof of well posedness for the original Cauchy problem. The multidimensional case is thus completely different.

3. What our example points out is that, to get well posedness, the initial data should be chosen in a space of functions smaller than $`𝐋^{\mathrm{}}`$. In view of the one-dimensional theory, it is natural to conjecture that a global weak solution should exist for initial data $`\overline{u}`$ with bounded variation. A proof of this conjecture (far from easy!) could rely on a compactness property for fluxes generated by O.D.E’s with right hand side in BV. More precisely, consider a sequence of smooth maps $`f_\nu :[0,T]\times IR^mIR^m`$ such that

$$\left|f_\nu (t,x)\right|C_1,$$
$`(3.4)`$
$$f_\nu _{BV}_0^T_{IR^m}\left|\frac{}{t}f_\nu \right|+\underset{i=1}{\overset{m}{}}\left|\frac{}{x_i}f_\nu \right|dxdtC_2.$$
$`(3.5)`$
Call $`tx(t)\mathrm{\Phi }_t^\nu (y)`$ the solution of

$$\dot{x}=f_\nu (t,x),x(0)=y.$$

Moreover, assume that the fluxes $`\mathrm{\Phi }_t^\nu `$ are all nearly incompressible, so that, for every bounded set $`AIR^m`$,

$$\frac{1}{C_3}\text{meas}(A)\text{meas}\left(\mathrm{\Phi }_t^\nu (A)\right)C_3\text{meas}(A),$$
$`(3.6)`$
for some constant $`C_3`$ and all $`t[0,T]`$, $`\nu 1`$. Conjecture: By possibly exstracting a subsequence, one has the convergence

$$\mathrm{\Phi }^\nu \mathrm{\Phi }\text{in}𝐋_{loc}^1$$

for some measurable flux $`\mathrm{\Phi }`$, also satisfying (3.6).

Notice that, by a well known compactness theorem, for some $`f𝐋^{\mathrm{}}`$ we certainly have

$$f_\nu f\text{in}𝐋_{loc}^1.$$

The existence of a limit flux $`\mathrm{\Phi }`$ would provide a new existence result concerning the Cauchy problem for the discontinuous O.D.E.

$$\dot{x}=f(t,x),x(0)=y$$
$`(3.7)`$
valid for a.e. initial data $`yIR^m`$, somewhat extending the result in .

References

 A.Bressan, Hyperbolic Systems of Conservation Laws. The One Dimensional Cauchy Problem, Oxford University Press, 2000.

 C. Dafermos, Hyperbolic Conservation Laws in Continuum Physics, Springer-Verlag, Berlin 1999.

 R. DiPerna and P. L. Lions, Ordinary differential equations, transport theory and Sobolev spaces, Invent. Math. 98 (1989), 511-517.

 S. Kruzhkov, First-order quasilinear equations with several space variables, Math. USSR Sbornik 10 (1970), 217–273.

 D. Serre, Systems of Conservation Laws I, II, Cambridge University Press, 2000.


text:
# 𝐵_{𝑠,𝑑}→𝛾𝛾 Decay with the Fourth Generation

## I Introduction

The minimal Standard Model (SM) is a very successful theory of the elementary particles known today. But it should not be the final theory, because it has too many unknown parameters to be put by hand. Most of these parameters are in the fermion part of the theory. We don’t know the origin of the quarks and leptons, as well as how to determinate their mass and number theoretically. We have to get their information all from experiments. There is still not a successful theory which can unify all the four basic interactions. Perhaps our elementary particles have substructure and we need a more fundamental theory than SM. This is beyond the scope of our current experiments. On the other hand, the recent measurement of the muon anomalous magnetic moment by the experiment E821  disagrees with the SM expectations at more than 1.6$`\sigma `$ \- 3$`\sigma `$  level. Moreover, there are convincing evidences that neutrinos are massive and oscillate in flavor . It seems to indicate the presence of new physics. The discovery of neutrino oscillation has been one of the most exciting experimental results in the recent years . As the experimental data accumulating in the Super-Kamiokande collaboration , more results published by SNO , K2K  and CHOOZ  experiments. These experiments give possible signals of new physics beyond SM.

In this paper, we consider the sequential fourth generation model to estimate the possible contributions to the decay $`B_{s,d}\gamma \gamma `$ from the exchange of the fourth up-type quark $`t^{}`$. From the point of phenomenology, there is a realistic question one may ask: What is the number of the fermion generations or weather there are additional quarks or leptons other than the known three generations. The present experiments can tell us there are only three generation of fermions with $`light`$ neutrinos whose mass is smaller than $`M_Z/2`$ , but the experiments don’t exclude the existence of other kind of additional generation, such as the fourth generation, with a $`heavy`$ neutrino, i.e. $`m_{\nu _4}M_Z/2`$ . Latest studies in the electroweak sector allow the existence of a fourth generation with heavy Dirac neutrino . Many authors studied models which extend the fermions part, such as vector-like quark models , sterile neutrino models  and the sequential four generation standard model (SM4) . Here we consider the sequential fourth generation, non -supersymmetric model, which includes an up-like quark $`t^{}`$, a down-like quark $`b^{}`$, a lepton $`\tau ^{}`$, and a heavy neutrino $`\nu ^{}`$ besides the known three generation of fermions in the SM. The properties of these new fermions are the same as their corresponding counterparts of other three generations except their masses and CKM matrix elements.

Many attempts have been made recently about the fourth generation  and the experimental searches of the fourth generation particles . In our previous papers , we investigated the rare $`B`$ meson decays with the fourth generation , $`ϵ^{^{}}/ϵ`$ in $`K^0`$ systems and $`\mathrm{\Delta }M_{B_{s,d}}`$ in SM4 . We got some interesting results, such as the new effects of the 4th generation particle on the meson decays and CP violation. We also got the constraints of the fourth generation CKM matrix factors . Recently, these decays in fourth generation were considered again. In ref. , the phases of the fourth generation CKM matrix factors, $`V_{t^{^{}}s}^{}V_{t^{^{}}b}`$, have been included. The exclusive B decays with the fourth generation were investigated in Ref. . These rare decays are very interesting to study the new physics effects of the fourth generation. In this paper, we would like to investigate the rare decay $`B_{s,d}\gamma \gamma `$ with the fourth generation.

As a loop-induced flavor changing neutral current (FCNC) process the inclusive decay (at quark level) $`bs\gamma \gamma `$ is in particular sensitive to contributions from those new physics beyond the SM. There is a vast interest in this decay. On the experimental side, only upper limits ($`90\%C.L.`$) on the branching ratios of $`B_{s,d}\gamma \gamma `$ decay are currently available

$`\mathrm{BR}(B_s\gamma \gamma )`$ $`<`$ $`1.48\times 10^4,\text{[25]},`$ (1)
$`\mathrm{BR}(B_d\gamma \gamma )`$ $`<`$ $`1.7\times 10^6,\text{[26]},`$ (2)
which are roughly two orders above the SM predictions ,

$`5\times 10^7<\mathrm{BR}(B_s\gamma \gamma )<15\times 10^7,`$ (3)
$`1\times 10^8<\mathrm{BR}(B_d\gamma \gamma )<4.0\times 10^8.`$ (4)
Within the SM, the electroweak contributions to $`bs\gamma \gamma `$ and $`B\gamma \gamma `$ decays have been calculated some time ago , the leading-order QCD corrections and the long-distance contributions were evaluated recently by several groups . The new physics corrections were also considered, for example, in the two-Higgs doublet model  and the supersymmetric model . There is also study concerning other exclusive channels . In a previous paper , we calculated the Technicolor corrections to the rare decays $`B_{s,d}\gamma \gamma `$ in the QCD factorization approach. The great progress in theoretical studies recently encourage us to do more investigations about this decay in the sequential fourth generation model.

This paper is organized as the following: In Section 2, we present the formulae for the decays $`B_{s,d}\gamma \gamma `$ in SM and the fourth generation model. We also obtain the constraints of the fourth generation CKM matrix factors, $`V_{t^{}d}^{}V_{t^{}b}`$ and $`V_{t^{}s}^{}V_{t^{}b}`$ from experiments in this section. In Section 3, we give the numerical analysis of the branching ratios in the fourth generation model. The conclusions are in the last section.

## II  $`B_{s,d}\gamma \gamma `$ in the fourth generation model

Same as in our last paper , up to the corrections of order $`1/m_W^2`$, the effective Hamiltonian for $`bs\gamma \gamma `$ at scales $`\mu _b=𝒪(m_b)`$ is just the one for $`bs\gamma `$ and takes the form

$$_{\mathrm{eff}}=\frac{G_\mathrm{F}}{\sqrt{2}}V_{ts}^{}V_{tb}\left[\underset{i=1}{\overset{6}{}}C_i(\mu _b)Q_i+C_{7\gamma }(\mu _b)Q_{7\gamma }+C_{8G}(\mu _b)Q_{8G}\right],$$
(5)
where in view of $`|V_{us}^{}V_{ub}/V_{ts}^{}V_{tb}|<0.02`$, we have neglected the term proportional to $`V_{us}^{}V_{ub}`$. Here $`Q_1\mathrm{}Q_6`$ are the usual four-fermion operators whose explicit form is given below. The last two operators in the Eq.(5), characteristic for this decay, are the magnetic–penguin operators. The complete list of operators is given as follows

$`O_1`$ $`=`$ $`(\overline{c}_{L\beta }\gamma ^\mu b_{L\alpha })(\overline{s}_{L\alpha }\gamma _\mu c_{L\beta }),`$ (6)
$`O_2`$ $`=`$ $`(\overline{c}_{L\alpha }\gamma ^\mu b_{L\alpha })(\overline{s}_{L\beta }\gamma _\mu c_{L\beta }),`$ (7)
$`O_3`$ $`=`$ $`(\overline{s}_{L\alpha }\gamma ^\mu b_{L\alpha }){\displaystyle \underset{q=u,d,s,c,b}{}}(\overline{q}_{L\beta }\gamma _\mu q_{L\beta }),`$ (8)
$`O_4`$ $`=`$ $`(\overline{s}_{L\alpha }\gamma ^\mu b_{L\beta }){\displaystyle \underset{q=u,d,s,c,b}{}}(\overline{q}_{L\beta }\gamma _\mu q_{L\alpha }),`$ (9)
$`O_5`$ $`=`$ $`(\overline{s}_{L\alpha }\gamma ^\mu b_{L\alpha }){\displaystyle \underset{q=u,d,s,c,b}{}}(\overline{q}_{R\beta }\gamma _\mu q_{R\beta }),`$ (10)
$`O_6`$ $`=`$ $`(\overline{s}_{L\alpha }\gamma ^\mu b_{L\beta }){\displaystyle \underset{q=u,d,s,c,b}{}}(\overline{q}_{R\beta }\gamma _\mu q_{R\alpha }),`$ (11)
$`O_7`$ $`=`$ $`(e/16\pi ^2)m_b\overline{s}_L\sigma ^{\mu \nu }b_RF_{\mu \nu },`$ (12)
$`O_8`$ $`=`$ $`(g/16\pi ^2)m_b\overline{s}_L\sigma ^{\mu \nu }T^ab_RG_{\mu \nu }^a.`$ (13)
It is the magnetic $`\gamma `$-penguin operator $`O_7`$, which plays the crucial role in this decay. The effective Hamiltonian for $`bd\gamma \gamma `$ is obtained from Eqs.(5-13) by the replacement $`sd`$.

In the case of the fourth generation there is an additional contribution to $`b(s,d)\gamma \gamma `$ from the virtual exchange of the fourth generation up quark $`t^{^{}}`$. In Fig.1, we draw the relevant Feynman diagrams which contribute to the decays $`b(s,d)\gamma \gamma `$. The Wilson coefficients of the dipole operators are given by

$$C_{7,8}^{\mathrm{eff}}(\mu _b)=C_{7,8}^{(\mathrm{SM})\mathrm{eff}}(\mu _b)+\frac{V_{t^{^{}}q}^{}V_{t^{^{}}b}}{V_{tq}^{}V_{tb}}C_{7,8}^{(4)\mathrm{eff}}(\mu _b),$$
(14)
where $`C_{7,8}^{(4)\mathrm{eff}}(\mu _b)`$ represent the contributions of $`t^{^{}}`$ to the Wilson coefficients. We recall here that the CKM coefficient corresponding to the $`t`$ quark contribution, i.e., $`V_{ts}^{}V_{tb}`$, is factorized in the effective Hamiltonian as given in Eq. (5). The formulae for calculating the Wilson coefficients $`C_{7,8}^{(4)}(m_W)`$ are the same as in the SM except exchanging $`t^{^{}}`$ quark for $`t`$ quark.

In SM4, the quark mixing matrix can be written as,

$$V=\left(\begin{array}{cccc}V_{ud}\hfill & V_{us}& \hfill V_{ub}& \hfill V_{ub^{}}\\ V_{cd}\hfill & V_{cs}& \hfill V_{cb}& \hfill V_{cb^{}}\\ V_{td}\hfill & V_{ts}& \hfill V_{tb}& \hfill V_{tb^{}}\\ V_{t^{}d}\hfill & V_{t^{}s}& \hfill V_{t^{}b}& \hfill V_{t^{}b^{}}\end{array}\right),$$
(15)
where $`V_{qb^{}}`$ and $`V_{t^{}q}`$ are the elements of the $`4\times 4`$ CKM mixing matrix of the SM4, which now contains nine parameters, i.e., six angles and three phases. The rest elements without prime in the matrix are the usual three generation CKM matrix elements.

The current experimental bounds of $`V_{t^{}s}^{}V_{t^{}b}`$ and $`V_{t^{}d}^{}V_{t^{}b}`$ can be found in ref. . The factor $`V_{t^{}s}^{}V_{t^{}b}`$ can be constrained by the decay $`BX_s\gamma `$

$`11.6\times 10^2<V_{t^{}s}^{}V_{t^{}b}^{(1)}<6.1\times 10^2,`$ (16)
or

$`1.9\times 10^3<V_{t^{^{}}s}^{}V_{t^{^{}}b}^{(2)}<3.6\times 10^3.`$ (17)
Considering the unitarity conditions of the CKM matrix,

$`{\displaystyle \underset{i}{}}V_{is}^{}V_{ib}=0,(i=u,c,t,t^{}),`$ (18)
and taking the average values of the SM $`3\times 3`$ CKM matrix elements from Ref. , we obtain,

$`|V_{t^{}s}^{}V_{t^{}b}|<7.6\times 10^2.`$ (19)
Combining the above constraint, we arrive at

$`7.6\times 10^2<V_{t^{}s}^{}V_{t^{}b}^{(1)}<6.0\times 10^2,`$ (20)
$`1.9\times 10^3<V_{t^{}s}^{}V_{t^{}b}^{(2)}<3.6\times 10^3.`$ (21)
In the following numerical calculation, we take the values $`V_{t^{}s}^{}V_{t^{}b}^{(1)}=7.0\times 10^2`$ and $`V_{t^{}s}^{}V_{t^{}b}^{(2)}=2.5\times 10^3`$, for illustration.

For the CKM factor $`V_{t^{}d}^{}V_{t^{}b}`$, we can get its constraint from the present experimental value of $`\mathrm{\Delta }M_{B_d}`$ . Also considering the unitarity conditions of the CKM matrix,

$`{\displaystyle \underset{i}{}}V_{id}^{}V_{ib}=0,(i=u,c,t,t^{}),`$ (22)
we finally get

$`1.0\times 10^4<V_{t^{}d}^{}V_{t^{}b}<0.5\times 10^4.`$ (23)
In the next section, we will take the values $`V_{t^{}d}^{}V_{t^{}b}=1.0\times 10^4`$, $`0.5\times 10^4`$, $`0.1\times 10^4`$, $`0.1\times 10^4`$ and $`0.5\times 10^4`$ in numerical calculation respectively.

## III The $`B_{s,d}\gamma \gamma `$ decay rates and phenomenology

We first calculate the decay $`B_s\gamma \gamma `$ (the formulas in the case of $`B_d\gamma \gamma `$ can be obtained just by the replacement of $`sd`$). For this exclusive decay of $`B_s`$ meson, we need to deal with both of the long distance and short distance QCD corrections in the $`B_s`$ meson side. Applying the factorization assumption, we shall adapt a phenomenological approach where the long distance effects are replaced by a few non-perturbative parameters. In other words we simply evaluate the hadronic matrix element of $`M_I+M_R`$ from one particle irreducible (1PI) and one particle reducible (1PR) Feynman diagrams, relying on a phenomenological model. Before doing so, it is important to note that $`M_R`$ is apparently non-local due to internal $`b`$ or $`s`$ quark propagator. To handle these non-local terms, one observes that the $`b`$ quark inside the $`B_s`$ meson carries most of the meson energy, and its four velocity can be treated as equal to that of $`B_s`$. Hence one may write $`b`$ quark momentum as $`p=m_bv`$ where $`v`$ is the common four velocity of $`b`$ and $`B_s`$. With this parametrization, we have

$`pk_1=m_bvk_1={\displaystyle \frac{1}{2}}m_bm_{B_s}=pk_2,`$ (24)
$`p^{}k_1=(pk_1k_2)k_1={\displaystyle \frac{1}{2}}m_{B_s}(m_{B_s}m_b)=p^{}k_2`$ , (25)
where the second equation is based on a constituent picture  that $`b`$ and $`\overline{s}`$ quarks share the total energy of $`B_s`$***Note that the momentum of $`\overline{s}`$ quark is $`p^{}`$ as $`p^{}`$ denotes the momentum of $`s`$ quark in $`bs\gamma \gamma `$.. Therefore $`p^{}`$ should be taken as the four momentum of the constituent $`\overline{s}`$ quark. With Eq.(25), $`M_R`$ is readily made local. We then compute the amplitude for $`B_s\gamma \gamma `$ using the following relations

$`0|\overline{s}\gamma _\mu \gamma _5b|B_s(P)`$ $`=`$ $`if_{B_s}P_\mu ,`$ (26)
$`0|\overline{s}\gamma _5b|B_s(P)`$ $`=`$ $`if_{B_s}M_B,`$ (27)
where $`f_{B_s}`$ is the $`B_s`$ meson decay constant which is about $`200`$ MeV according to recent Lattice QCD calculations .

The total amplitude is now separated into a CP-even and a CP-odd part

$$T(B_s\gamma \gamma )=M^+F_{\mu \nu }F^{\mu \nu }+iM^{}F_{\mu \nu }\stackrel{~}{F}^{\mu \nu }.$$
(28)
We find that

$$M^+=\frac{4\sqrt{2}\alpha G_F}{9\pi }f_{B_s}V_{ts}^{}V_{tb}\left(\frac{m_b}{m_{B_s}}BK(m_b^2)+\frac{3C_7}{8\overline{\mathrm{\Lambda }}}\right),$$
(29)
with $`B=(3C_6+C_5)/4`$, $`\overline{\mathrm{\Lambda }}=1m_b/m_{B_s}`$, and

$$M^{}=\frac{4\sqrt{2}\alpha G_F}{9\pi }f_{B_s}V_{ts}^{}V_{tb}\left(\underset{q}{}A_qJ(m_q^2)+\frac{m_b}{m_{B_s}}BL(m_b^2)+\frac{3C_7}{8\overline{\mathrm{\Lambda }}}\right),$$
(30)
where

$`A_u`$ $`=`$ $`(C_3C_5)N_c+(C_4C_6),`$ (31)
$`A_d`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[(C_3C_5)N_c+(C_4C_6)\right],`$ (32)
$`A_c`$ $`=`$ $`(C_1+C_3C_5)N_c+(C_2+C_4C_6),`$ (33)
$`A_s`$ $`=`$ $`A_b={\displaystyle \frac{1}{4}}\left[(C_3+C_4C_5)N_c+(C_3+C_4C_6)\right].`$ (34)
The functions $`J(m^2)`$, $`K(m^2)`$ and $`L(m^2)`$ are defined by

$`J(m^2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{z}}I_{00}(m^2),`$ (35)
$`K(m^2)`$ $`=`$ $`2+{\displaystyle \frac{4z}{z}}I_{00}(m^2),`$ (36)
$`L(m^2)`$ $`=`$ $`I_{00}(m^2),`$ (37)
where $`z=m_{B_s}^2/m^2`$, and

$`I_{00}(m^2)`$ $`=`$ $`{\displaystyle _0^1}𝑑x{\displaystyle \frac{1}{x}}\mathrm{ln}(1zx+zx^2)`$ (38)
$`=`$ $`\{\begin{array}{cc}2\mathrm{arctan}^2(\sqrt{\frac{z}{4z}});\hfill & z<4\hfill \\ \frac{\pi ^2}{2}2\mathrm{log}^2(\frac{\sqrt{z}+\sqrt{z4}}{2})+2i\pi \mathrm{log}(\frac{\sqrt{z}+\sqrt{z4}}{2});\hfill & z>4\hfill \end{array}.`$ (41)
The decay width for $`B_s\gamma \gamma `$ is simply

$$\mathrm{\Gamma }(B_s\gamma \gamma )=\frac{m_{B_s}^3}{16\pi }(|M^+|^2+|M^{}|^2).$$
(42)
To obtain numerical results, we have set light quark masses to zero and used  $`m_t=175\mathrm{GeV}`$, $`m_b=5.0\mathrm{GeV}`$ and $`m_c=1.5\mathrm{GeV}`$. Furthermore, we take $`m_{B_s}=5.37\mathrm{GeV}`$, $`m_{B_d}=5.28\mathrm{GeV}`$ and $`\alpha =1/129`$. The numerical values for Wilson coefficients $`C_1C_8`$ evaluated in SM at $`\mu =m_b`$ are listed in Table I.

We first consider the decay mode $`B_s\gamma \gamma `$. Figs.2, 3 and 4 show the branching ratio $`Br(B_s\gamma \gamma )`$ $`versus`$ the mass of $`t^{}`$ and the fourth generation CKM matrix factors $`V_{t^{}s}^{}V_{t^{}b}^{(1)}`$ and $`V_{t^{}s}^{}V_{t^{}b}^{(2)}`$, respectively. From these figures, we find that the SM4 predictions for the branching ratio $`Br(B_s\gamma \gamma )`$ are rather different from that of the standard model. The decay rates can be enhanced by about 1-2 orders relative to the SM prediction in the reasonable range of $`m_t^{}`$ and $`V_{t^{}s}^{}V_{t^{}b}`$. This implies that the new physics signal of the fourth generation Model may come from the decay of $`B_s\gamma \gamma `$. The future experiment of LHCb in CERN can test it. From Fig.2, we find that the branching ratio of $`B_s\gamma \gamma `$ increase with the mass of $`t^{}`$ for the fourth generation CKM matrix factor, $`V_{t^{}s}^{}V_{t^{}b}=7.0\times 10^2`$ and $`2.5\times 10^3`$. When $`V_{t^{}s}^{}V_{t^{}b}=7.0\times 10^2`$, the curve increase rapidly with $`m_t^{}`$. But when $`V_{t^{}s}^{}V_{t^{}b}=2.5\times 10^3`$, the branching ratio changes very slowly. Because the fourth generation CKM matrix factor, $`V_{t^{}s}^{}V_{t^{}b}`$, is small in the latter case. It can be seen from Eq. (14) that the new physics contribution is proportional to the size of the CKM factor $`V_{t^{}s}^{}V_{t^{}b}`$, thus it is very important in the new physics contribution with $`m_t^{}`$. From Figs.3 and 4, one can also see that, the branching ratio of $`B_s\gamma \gamma `$ decease with $`V_{t^{}s}^{}V_{t^{}b}^{(1)}`$ but increase with $`V_{t^{}s}^{}V_{t^{}b}^{(2)}`$, since $`V_{t^{}s}^{}V_{t^{}b}^{(1)}`$ is negative and gives the opposite contributions as SM to the coefficients $`C_{7,8}^{eff}`$ (see Eq. (14) ).

For the decay $`B_d\gamma \gamma `$, it is easy to find from Figs.5 and 6 that the new physics contribution to the branching ratio is tiny. It is shown that it is not suitable to test the new physics signals from the 4th generation in this channel. The prediction of the branching ratio in the 4th generation model is between $`2.35\times 10^8`$ and $`2.45\times 10^8`$, when we take the reasonable values of $`m_t^{}`$ (from 200 GeV to 1000 GeV) and $`V_{t^{}d}^{}V_{t^{}b}`$ (see Eq. (23)). This numerical result is fully in the range of the SM prediction (see Eq. (4)). The reason is that the experimental constraint of the 4th generation CKM matrix element factor $`V_{t^{}d}^{}V_{t^{}b}`$ is very stringent (only at $`10^4`$ level, see Eq.(23)). We also find that the curves change little along with the variations of $`m_t^{}`$ and $`V_{t^{}d}^{}V_{t^{}b}`$. The reason is also from the tiny $`V_{t^{}d}^{}V_{t^{}b}`$. We can see that the branching ratios of both $`B_s\gamma \gamma `$ and $`B_s\gamma \gamma `$ decays strongly depend on the size of the related mixing elements.

## IV Conclusion

As a summary, the size of new physics contributions to the rare decays of $`B_{s,d}\gamma \gamma `$ from the fourth generation strongly depends on the values of the masses of the new heavy up quark $`t^{}`$ and the 4th generation CKM matrix elements, $`V_{t^{}s}^{}V_{t^{}b}`$ and$`V_{t^{}d}^{}V_{t^{}b}`$. The results in the 4th generation model are quite different from the SM case for decay $`B_s\gamma \gamma `$, its branching ratio can be enhanced by about 1-2 orders of magnitude. A large branching ratio of $`B_s\gamma \gamma `$ decay measured in experiment may be a new physics signal of the 4th generation model. But for decay $`B_d\gamma \gamma `$, because of the tiny values of 4th generation CKM matrix element, $`V_{t^{}d}^{}V_{t^{}b}`$, the new physics contribution to $`\mathrm{BR}(B_d\gamma \gamma )`$ is very small and can not be distinct from the SM prediction.

ACKNOWLEDGMENT

W.J. Huo acknowledges the support from the Chinese Postdoctoral Science Foundation.from the Natural Sciences and Engineering Research Council of Canada. C.D. Lü acknowledges the support by National Science Foundation of China under Grants No. 90103013 and 10135060. And Z.J. Xiao acknowledges the support by the National Natural Science Foundation of China under Grants No. 10075013 and 10275035, and by the Research Foundation of Nanjing Normal University under Grant No. 214080A916.


text:
# Forward-backward Asymmetry in 𝐾⁺→𝜋⁺ℓ⁺ℓ⁻

## I Introduction

The flavor-changing neutral current (FCNC) processes of $`K^\pm \pi ^\pm \mathrm{}^+\mathrm{}^{}(\mathrm{}=e,\mu )`$ are suppressed and dominated by the long distance (LD) contributions involving one photon exchange k1 ; k2 ; k3 ; k4  in the standard model (SM). The decays have been successfully described within the framework of chiral perturbation theory (ChPT) ChPT0  including electroweak interactions at $`𝒪(p^6)`$ ChPT1  in terms of a vector interaction form factor fixed by experiments. However, it is important to compare the measurements in the two decays to see if there are differences in the form factors since they would indicate new physics. Recently, the vector form factor has been determined by the high precision measurement on the electron mode by the BNL-E865 Collaboration E865e  at the Brookhaven Alternating Gradient Synchrotron (AGS) with a sample of 10300 events and branching ratio (BR) of $`[2.94\pm 0.05(stat)\pm 0.13(syst)\pm 0.05(theor)]\times 10^7`$. For the muon channel, it was first observed by BNL-E787 E787  at AGS with the measured branching ratio being $`(5.0\pm 0.4\pm 0.9)\times 10^8`$, which is too small to accommodate within the SM. However, two subsequent experiments of BNL-E865 E865m  and HyperCP (E871) E871  have measured that $`BR(K^+\pi ^+\mu ^+\mu ^{})=(9.22\pm 0.60\pm 0.49)`$ and $`(9.8\pm 1.0\pm 0.5)\times 10^8`$, respectively, which are all consistent with the model-independent analysis based on the data of $`BR(K^+\pi ^+e^+e^{})`$. However, there are still rooms for new physics, particularly in the muon mode.

In Refs. geng1 ; geng2 , $`P`$ and $`T`$ violating muon polarization effects in $`K^+\pi ^+\mu ^+\mu ^{}`$ were discussed in various theoretical models. In this paper, we study the forward-backward asymmetry (FBA) in the decay of $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$ with $`\mathrm{}=e`$ or $`\mu `$. It is known that the FBA in $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$ violates $`P`$ like the longitudinal lepton polarization but it vanishes in the SM and can only exist if there is a scalar type interaction. In the multi-Higgs doublet models such as the most popular two-Higgs doublet (2HDM) of type II 2HDM , where two Higgs scalar doublets ($`H_u`$ and $`H_d`$) are coupled to up- and down-type quarks, respectively, the scalar type of four fermion operators $`\overline{s}_Rd_L\overline{\mathrm{}}\mathrm{}`$ can be generated at the loop level tanB0 ; tanB . This type of operators is particularly interesting in the minimal supersymmetric standard model (MSSM) MSSM  since it receives an enormous enhancement for the large ratio of $`v_u/v_d=\mathrm{tan}\beta `$ where $`v_{u(d)}`$ is the vacuum expectation value of the Higgs doublet $`H_{u(d)}`$. Recently, there has been considerable interest for the large $`\mathrm{tan}\beta `$ effects in $`B`$ decays such as $`B\mu ^+\mu ^{}`$ and $`BK\mu ^+\mu ^{}`$ tanB0 ; tanB ; Bmm . In the report, we will discuss the large $`\mathrm{tan}\beta `$ scenario in the MSSM for $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$.

The paper is organized as follows. In Sec. II, we present the general analysis for the forward-backward asymmetries in $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$ ($`\mathrm{}=e`$ and $`\mu `$). In Sec. III, we discuss the experimental constraints on the asymmetries. We estimate the asymmetries in the MSSM in Sec. IV. We present our conclusions in Sec. V.

## II General analysis

We write the decay as

$`K^+(p_K)\pi ^+(p_\pi )\mathrm{}^+(p)\mathrm{}^{}(\overline{p}),`$ (1)
where $`p_K`$, $`p_\pi `$, $`p`$ and $`\overline{p}`$ are four-momenta of $`K^+`$, $`\pi ^+`$, $`\mathrm{}^+`$ and $`\mathrm{}^{}`$, respectively. The most general invariant amplitude for the decay can be written as geng1 ; geng2 ; geng3

$``$ $`=`$ $`F_S\overline{\mathrm{}}\mathrm{}+iF_P\overline{\mathrm{}}\gamma _5\mathrm{}+F_Vp_K^\mu \overline{\mathrm{}}\gamma _\mu \mathrm{}+F_Ap_K^\mu \overline{\mathrm{}}\gamma _\mu \gamma _5\mathrm{},`$ (2)
where $`F_S`$, $`F_P`$, $`F_V`$ and $`F_A`$ are scalar, pseudo-scalar, vector and axial-vector form factors, respectively. The differential decay rate in the $`K^+`$ rest frame is given by geng1

$`{\displaystyle \frac{d^2\mathrm{\Gamma }}{dEd\overline{E}}}`$ $`=`$ $`{\displaystyle \frac{1}{2^4\pi ^3m_K}}\left[\right|F_S|^2{\displaystyle \frac{1}{2}}(s4m_l^2)+|F_P|^2{\displaystyle \frac{1}{2}}s`$ (3)
$`+\left|F_V\right|^2m_K^2(2E\overline{E}{\displaystyle \frac{1}{2}}s)+\left|F_A\right|^2m_K^2(2E\overline{E}{\displaystyle \frac{1}{2}}s+2m_l^2)`$
$`+2\text{Re}(F_SF_V^{})m_lm_K(E\overline{E})+\text{Im}(F_PF_A^{})m_l(m_\pi ^2m_K^2s)],`$
where $`m_l`$ is the lepton mass, $`E(\overline{E})`$ is the energy of $`\mu ^+(\mu ^{})`$ and $`s=(p+\overline{p})^2=2(m_l^2+E\overline{E}𝐩\overline{𝐩})`$ is the invariant mass of the dilepton system. In terms of the invariant mass and the angle $`\theta `$ between the three-momentum of the kaon and the three-momentum of the $`\mathrm{}^{}`$ in the dilepton rest frame, we can rewrite Eq. (3) as

$`{\displaystyle \frac{d^2\mathrm{\Gamma }}{dsd\mathrm{cos}\theta }}`$ $`=`$ $`{\displaystyle \frac{1}{2^8\pi ^3m_K^3}}\beta _l\lambda ^{\frac{1}{2}}(s)\left\{\right|F_S|^2s\beta _l^2+|F_P|^2s`$ (4)
$`+\left|F_V\right|^2{\displaystyle \frac{1}{4}}\lambda (s)(1\beta _l^2\mathrm{cos}^2\theta )+\left|F_A\right|^2\left[{\displaystyle \frac{1}{4}}\lambda (s)(1\beta _l^2\mathrm{cos}^2\theta )+4m_K^2m_l^2\right]`$
$`+\text{Re}(F_SF_V^{})2m_l\beta _l\lambda ^{\frac{1}{2}}(s)\mathrm{cos}\theta +\text{Im}(F_PF_A^{})2m_l(m_\pi ^2m_K^2s)\},`$
where $`\lambda (s)=m_K^4+m_\pi ^4+s^22m_\pi ^2s2m_K^2s2m_\pi ^2m_K^2`$ and $`\beta _l=(14m_l^2/s)^{\frac{1}{2}}`$ with $`s`$ and $`\mathrm{cos}\theta `$ bounded by

$`4m_l^2s(m_Km_\pi )^2,1\mathrm{cos}\theta 1.`$ (5)
Here, we have used that

$`E={\displaystyle \frac{s+m_K^2m_\pi ^2+\beta _l\lambda ^{\frac{1}{2}}(s)\mathrm{cos}\theta }{4m_K}},\overline{E}={\displaystyle \frac{s+m_K^2m_\pi ^2\beta _l\lambda ^{\frac{1}{2}}(s)\mathrm{cos}\theta }{4m_K}}.`$ (6)
By integrating the angle $`\theta `$ in Eq. (4), we obtain

$`{\displaystyle \frac{d\mathrm{\Gamma }}{ds}}`$ $`=`$ $`{\displaystyle \frac{1}{2^8\pi ^3m_K^3}}\beta _l\lambda ^{\frac{1}{2}}(s)\left\{\right|F_S|^22s\beta _l^2+|F_p|^22s+|F_V|^2{\displaystyle \frac{1}{3}}\lambda (s)(1+{\displaystyle \frac{2m_l^2}{s}})`$ (7)
$`+|F_A|^2[{\displaystyle \frac{1}{3}}\lambda (s)(1+{\displaystyle \frac{2m_l^2}{s}})+8m_K^2m_l^2]+\text{Im}(F_PF_A^{})4m_l(m_\pi ^2m_K^2s)\}.`$
From Eq. (4) and the definition of the forward-backward asymmetry

$`A_{FB}(s)`$ $``$ $`{\displaystyle \frac{_0^1d\mathrm{cos}\theta \frac{d^2\mathrm{\Gamma }}{dsd\mathrm{cos}\theta }_1^0d\mathrm{cos}\theta \frac{d^2\mathrm{\Gamma }}{dsd\mathrm{cos}\theta }}{_0^1d\mathrm{cos}\theta \frac{d^2\mathrm{\Gamma }}{dsd\mathrm{cos}\theta }+_1^0d\mathrm{cos}\theta \frac{d^2\mathrm{\Gamma }}{dsd\mathrm{cos}\theta }}},`$ (8)
we find that

$`A_{FB}(s)`$ $`=`$ $`{\displaystyle \frac{1}{2^8\pi ^3m_K^3}}2m_l\beta _l^2\lambda (s)\text{Re}(F_SF_V^{})\left({\displaystyle \frac{d\mathrm{\Gamma }}{ds}}\right)^1.`$ (9)
As seen from Eq. (9), to get a nonzero value of $`A_{FB}`$, it is necessary to have a scalar interaction. However, in the SM the contributions from $`F_S`$ to the decay widths of $`K^+\pi ^+e^+e^{}`$ and $`K^+\pi ^+\mu ^+\mu ^{}`$ are about 7 and 4 orders of magnitude smaller than those from $`F_V`$ ChPT2 ; Gao , respectively, and therefore the forward-backward asymmetries are expected to be vanishingly small.

## III Experimental constraints

To study the experimental constraints on $`A_{FB}`$ in $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$, we consider the amplitude adopted in Ref. E865e ;

$``$ $`=`$ $`{\displaystyle \frac{\alpha G_F}{4\pi }}f_VP^\mu \overline{\mathrm{}}\gamma _\mu \mathrm{}+G_Fm_Kf_S\overline{\mathrm{}}\mathrm{}+G_Ff_T{\displaystyle \frac{P^\mu q^\nu }{m_K}}\overline{\mathrm{}}\sigma _{\mu \nu }\mathrm{}`$ (10)
where $`f_{V,S,T}`$ are dimensionless form factors of vector, scalar, and tensor interactions, respectively, $`P=p_K+p_\pi `$ and $`q=p_Kp_\pi `$. It is clear that, in Eq. (10), the vector term arises from the one photon exchange in the SM, which gives the dominant contribution to the decay rate, whereas the scalar and tensor ones from some new physics beyond the SM Tensor .

For the form factor $`f_V`$, we take the form derived in the ChPT ChPT1 , given by

$`f_V(s)`$ $`=`$ $`a_++b_+{\displaystyle \frac{s}{m_K^2}}+\omega ^{\pi \pi }(s),`$ (11)
where $`a_+`$ and $`b_+`$ are free parameters and $`\omega ^{\pi \pi }`$ is the contribution from a pion loop diagram given in Ref. ChPT1 . The experimental measurement on $`K^+\pi ^+e^+e^{}`$ at BNL-E865 E865e  has determined the parameters of $`a_+`$ and $`b_+`$ to be $`0.587\pm 0.010`$ and $`0.655\pm 0.044`$, respectively. The scalar and tensor form factors in Eq. (10) for $`K^+\pi ^+e^+e^{}`$ are also constrained by the experiment E865e  and the results are that

$`|f_S|`$ $`<`$ $`6.6\times 10^5or|f_T|<\mathrm{\hspace{0.33em}3.7}\times 10^4`$ (12)
for the existence of either scalar or tensor interaction. We note that so far there are no similar constraints on $`f_{S,T}`$ for $`K^+\pi ^+\mu ^+\mu ^{}`$ and they can be quite different for the two channels in theoretical models.

It is easy to see that the amplitude in Eq. (10) can be simplified to

$``$ $`=`$ $`{\displaystyle \frac{\alpha G_F}{4\pi }}f_V^{}P^\mu \overline{\mathrm{}}\gamma _\mu \mathrm{}+G_Fm_Kf_s^{}\overline{\mathrm{}}\mathrm{}`$ (13)
with

$`f_V^{}`$ $`=`$ $`f_V{\displaystyle \frac{8\pi im_l}{\alpha m_K}}f_T,f_S^{}=f_S{\displaystyle \frac{i\beta _l\lambda ^{\frac{1}{2}}(s)\mathrm{cos}\theta }{m_K^2}}f_T.`$ (14)
By comparing the amplitude in Eq. (13) with the general one in Eq. (2), we get

$`F_V={\displaystyle \frac{\alpha G_F}{2\pi }}f_V^{},F_S=G_Fm_Kf_S^{},F_{P,A}=0.`$ (15)
From Eqs. (4), (7) and (15), we obtain

$`{\displaystyle \frac{d^2\mathrm{\Gamma }}{dsd\mathrm{cos}\theta }}`$ $`=`$ $`{\displaystyle \frac{G_F^2}{2^8\pi ^3m_K^3}}\beta _l\lambda ^{\frac{1}{2}}(s)\left\{\right|f_V|^2{\displaystyle \frac{\alpha ^2}{16\pi ^2}}\lambda (s)(1\beta _l^2\mathrm{cos}^2\theta )+|f_S|^2s\beta _l^2m_K^2`$ (16)
$`+|f_T|^2{\displaystyle \frac{s\lambda (s)}{m_K^2}}(\mathrm{cos}^2\theta +{\displaystyle \frac{4m_l^2}{s}}\mathrm{sin}^2\theta )+\text{Re}(f_V^{}f_S){\displaystyle \frac{\alpha m_lm_K}{\pi }}\beta _l\lambda ^{\frac{1}{2}}(s)\mathrm{cos}\theta `$
$`\text{Im}(f_Vf_T^{}){\displaystyle \frac{\alpha \lambda (s)}{\pi }}{\displaystyle \frac{m_l}{m_K}}\text{Im}(f_Sf_T^{})2s\beta \lambda ^{\frac{1}{2}}(s)\mathrm{cos}\theta \}`$
and

$`{\displaystyle \frac{d\mathrm{\Gamma }}{ds}}`$ $`=`$ $`{\displaystyle \frac{G_F^2}{2^8\pi ^3m_K^3}}\beta _l\lambda ^{\frac{1}{2}}(s)\left\{\right|f_V|^2{\displaystyle \frac{\alpha ^2\lambda (s)}{4\pi ^2}}{\displaystyle \frac{1}{3}}(1+{\displaystyle \frac{2m_l^2}{s}})+2|f_S|^2s\beta _l^2m_K^2`$ (17)
$`+|f_T|^2{\displaystyle \frac{2s\lambda (s)}{3m_K^2}}(1+{\displaystyle \frac{8m_l^2}{s}})\text{Im}(f_Vf_T^{}){\displaystyle \frac{2\alpha \lambda (s)}{\pi }}{\displaystyle \frac{m_l}{m_K}}\}.`$
Similarly, from Eq. (9) we find

$`A_{FB}(s)`$ $`=`$ $`{\displaystyle \frac{G_F^2}{2^8\pi ^3m_K^3}}\beta _l^2\lambda (s)\left[\text{Re}(f_V^{}f_S){\displaystyle \frac{\alpha m_lm_K}{\pi }}\text{Im}(f_Sf_T^{})2s\right]\left({\displaystyle \frac{d\mathrm{\Gamma }}{ds}}\right)^1.`$ (18)
From Eq. (17), one can check that the bound for $`f_S`$ or $`f_T`$ in Eq. (12) yields at most a few percent of the decay rate in $`K^+\pi ^+e^+e^{}`$. Moreover, the last term in Eq. (17) is negligible for the electron channel no matter whether $`f_T`$ is real or imaginary due to the electron mass suppression. However, for the muon case this term could be large and spoil the vector dominant mechanism if the imaginary part of $`f_T`$ is not small. In Figure 1, we show the differential decay rate and forward-backward asymmetry as functions of $`\widehat{s}=s/m_K^2`$ for the decay of $`K^+\pi ^+e^+e^{}`$ by using the upper value of $`f_S`$ in Eq. (12) and $`f_T=0`$. In Figure 2, we display them by assuming that $`f_s4\times 10^5i`$ and $`f_T2\times 10^4`$. As illustrations, in Figures 3 and 4 we also give $`d\mathrm{\Gamma }/d\widehat{s}`$ and $`A_{FB}`$ in $`K^+\pi ^+\mu ^+\mu ^{}`$ with the same sets of parameters as those in Figures 1 and 2, respectively. It is clear that, as mentioned early, since there is no direct strict experimental constraint on $`f_S`$ or $`f_T`$ in the muon mode, $`A_{FB}(K^+\pi ^+\mu ^+\mu ^{})`$ can be arbitrary large.

## IV Supersymmetry

In the MSSM, the one-loop effective down-type Yukawa interaction is given by

$`^{eff}`$ $`=`$ $`\overline{d}_RY_d[H_d+(ϵ_0+ϵ_YY_u^{}Y_u)H_u^{}]Q_L+h.c.,`$ (19)
where $`Y_{u,d}`$ are $`3\times 3`$ Yukawa coupling matrices and $`ϵ_{0,Y}`$ are defined in Ref. tanB , which are typically $`𝒪(10^2)`$. In the diagonal $`Y_d`$ basis of $`(Y_d)_{ij}=y_i^d\delta _{ij}`$, the interaction in Eq. (19) becomes

$`_{mass}^{eff}`$ $`=`$ $`v_d\overline{d}_R^iy_i^d\left[(1+ϵ_0\mathrm{tan}\beta )\delta _{ij}+ϵ_YV_{ik}^{}(y_k^u)^2V_{kj}\right]d_L^j+h.c.,`$ (20)
where $`V`$ is the CKM mixing matrix. By writing the effective Hamiltonian in the transition $`sd\mathrm{}^+\mathrm{}^{}`$ induced by the scalar type of interactions as

$`_S^{eff}`$ $`=`$ $`(C_S\overline{s}_Rd_L+C_S^{}\overline{d}_Rs_L)\overline{\mathrm{}}\mathrm{},`$ (21)
from Eq. (20) one has that tanB

$`C_S`$ $`=`$ $`{\displaystyle \frac{G_F^2}{4\pi ^2}}{\displaystyle \frac{m_sm_lm_t^2\overline{\lambda }_{21}^t\mathrm{tan}^3\beta }{(1+ϵ_0\mathrm{tan}\beta )}}{\displaystyle \frac{1}{M_A^2}}{\displaystyle \frac{\mu Af(x_{\mu L},x_{\mu R})}{M_{\stackrel{~}{t}_L}^2}},`$
$`C_S^{}`$ $``$ $`{\displaystyle \frac{m_d}{m_s}}C_S,`$ (22)
where $`A`$ is the coupling of the soft-breaking trilinear term and

$`\overline{\lambda }_{21}^t`$ $`=`$ $`\lambda _{21}^t\left[{\displaystyle \frac{1+\mathrm{tan}\beta (ϵ_0+ϵ_Yy_t^2)}{1+ϵ_0\mathrm{tan}\beta }}\right]^2`$
$`f(x,y)`$ $`=`$ $`{\displaystyle \frac{1}{xy}}\left[{\displaystyle \frac{x\mathrm{ln}x}{1x}}{\displaystyle \frac{y\mathrm{ln}y}{1y}}\right],f(1,1)={\displaystyle \frac{1}{2}},`$ (23)
with

$`x_{\mu L}`$ $`=`$ $`{\displaystyle \frac{\mu ^2}{M_{\stackrel{~}{t}_L}^2}},x_{\mu R}={\displaystyle \frac{M_{\stackrel{~}{t}_R}^2}{M_{\stackrel{~}{t}_L}^2}},\lambda _{21}^t=V_{ts}^{}V_{td},`$ (24)
and $`y_t`$ being the top quark Yukawa coupling. By comparing Eq. (21) with Eq. (10) and using

$`<\pi |\overline{d}(1+\gamma _5)s|K>`$ $``$ $`{\displaystyle \frac{m_K^2}{m_s}}f_+,`$ (25)
we find

$`f_S^{MSSM}`$ $``$ $`{\displaystyle \frac{G_F^2}{8\pi ^2}}{\displaystyle \frac{m_Km_lm_t^2\overline{\lambda }_{21}^t\mathrm{tan}^3\beta }{(1+ϵ_0\mathrm{tan}\beta )^2}}{\displaystyle \frac{1}{M_A^2}}{\displaystyle \frac{\mu Af(x_{\mu L},x_{\mu R})}{M_{\stackrel{~}{t}_L}^2}},`$ (26)
where we have neglected the small terms related to $`y_{1,2}^u`$ and used $`f_+1`$.

To estimate the scalar form factor in Eq. (26) in the MSSM with large $`\mathrm{tan}\beta `$, we take $`ϵ_01/100ϵ_Yy_t^2`$ and $`\mathrm{tan}\beta =50r`$ and we get

$`f_S^{MSSM}|_{\mathrm{}=e}`$ $``$ $`1.1\times 10^9(1\overline{\rho }i\overline{\eta }){\displaystyle \frac{r^3}{(1+\frac{1}{2}r)^2}}\left({\displaystyle \frac{200GeV}{M_A}}\right)^2\left({\displaystyle \frac{\mu Af(x_{\mu L},x_{\mu R})}{M_{\stackrel{~}{t}_L}^2}}\right),`$
$`f_S^{MSSM}|_{\mathrm{}=\mu }`$ $``$ $`2.3\times 10^7(1\overline{\rho }i\overline{\eta }){\displaystyle \frac{r^3}{(1+\frac{1}{2}r)^2}}\left({\displaystyle \frac{200GeV}{M_A}}\right)^2\left({\displaystyle \frac{\mu Af(x_{\mu L},x_{\mu R})}{M_{\stackrel{~}{t}_L}^2}}\right),`$ (27)
where $`\overline{\rho }=\rho (1\lambda ^2/2)`$ and $`\overline{\eta }=\eta (1\lambda ^2/2)`$ with $`\lambda `$, $`\rho `$ and $`\eta `$ being the Wolfenstein parameters of the CKM matrix $`V`$. Since the values of $`f_S^{MSSM}|_{\mathrm{}=e,\mu }`$ in Eq. (27) are about three and one orders of magnitude smaller than the experimental bound in Eq. (12) and thus the scalar contributions to the decay rates in the MSSM are negligible, respectively. Moreover, the scalar contribution to the FBA in $`K^+\pi ^+e^+e^{}`$ is also suppressed. However, in $`K^+\pi ^+\mu ^+\mu ^{}`$ the FBA can be as large as $`10^3`$ as shown in Figure 5 by using $`\overline{\rho }0.2`$ Buras02  and assuming $`r1`$, $`M_A200GeV`$ and $`\mu Af(x_{\mu L},x_{\mu R})/M_{\stackrel{~}{t}_L}^22`$. We note that $`A_{FB}(K^+\pi ^+\mu ^+\mu ^{})=𝒪(10^3)`$ is accessible to future experiments such as the CKM at Fermilab, where the order of $`10^5`$ events can be produced Kaon .

## V Conclusions

We have studied the forward-backward asymmetries in the decays of $`K^+\pi ^+\mathrm{}^+\mathrm{}^{}`$ ($`\mathrm{}=e`$ and $`\mu `$) in the most general amplitudes. In particular, we have explored the experimental constraints on the asymmetries by including the scalar and tensor interactions. We have found that with the scalar (tensor) term the asymmetry can be up to $`𝒪(10^3)`$ ($`𝒪(10^1))`$ and arbitrary large for the electron and muon channels, respectively, without conflict with the experimental data. We have also discussed the asymmetries in the minimal supersymmetric standard model where the scalar terms can be explicitly induced. We have shown that the FBA in $`K^+\pi ^+e^+e^{}`$ is negligibly small due to the electron mass suppression, but in $`K^+\pi ^+\mu ^+\mu ^{}`$ it can be as large as $`𝒪(10^3)`$ with the large $`\mathrm{tan}\beta `$, which can be tested in future experiments such as the CKM experiment at Fermilab.

Acknowledgments

We thank Drs. D.N. Gao and H.K. Park for useful communications and Prof. T. Han for encouragement and discussions. This work was supported in part by the National Science Council of the Republic of China under Contract Nos. NSC-91-2112-M-001-053 and NSC-91-2112-M-007-043.


text:
# The Multiplicity Function of Galaxies

## 1 Introduction

The multiplicity function (hereafter MF), in its differential form, is defined as the number of groups or clusters per area or volume unit and per richness unit.
The MF, which is the richness spectrum of galaxy aggregates, parametrises the observed clustering of galaxies and hence, together with the correlation and luminosity functions, is one of the fundamental cosmological observables. With respect to the complete description of clustering, the MF is complementary to the covariance function (which is related to the two-point correlation function), being related to the ratio of the amplitude of the higher-order to the two-point correlation functions (Gott & Turner, 1977, hereafter GT). Due to computational costs and errors, the measurement of correlation functions of order $`N`$ becomes unreliable for $`N>3`$, and the MF is therefore a crucial means of obtaining information on higher order clustering.
The Press–Schechter theory (Press & Schecter, 1974) states that the shape of the mass function (a power-law mass distribution with an exponential cutoff at the bright end) should provide important clues concerning the conditions at the epoch of recombination and does not depend on the cosmic density parameter $`\mathrm{\Omega }`$. The steepness of the initial density fluctuation spectrum constrains the broadness of the mass function.
The MF, the mass function or the luminosity function all describe in a similar way the cosmic abundance of objects and, in fact, present similar shape (Bachall, 1979).

Despite the fact that the early descriptions of galaxy clustering properties were given in terms of the MF (Gott & Turner, 1977), most authors have focused on the shape of the mass function, which can be directly compared to the PS formalism. Even when the observed quantity is the MF, some authors (Bachall & Cen, 1993) prefer to convert it into a mass function using a reliable M/L ratio. Nevertheless, one must consider all of the uncertainties introduced by the mass estimation, which are propagated to the mass function determination. These include errors in the internal velocity dispersion used for dynamical mass estimates, the large intrinsic scatter in the richness-mass relation, and errors in assuming dynamical equilibrium for all clusters when using X-ray data (Girardi et al., 1998).

The main problem which must be overcome when determining the MF is the production of a statistically significant and unbiased catalogue of groups and clusters covering a large enough area of the sky and encompassing cosmic structures spanning a wide range of richness, from very low multiplicity structures such as galaxy triplets, up to very rich clusters with several hundred members.

In the past, catalogues of groups and clusters have been derived from either 3-D data (cf. Maia et al., 1989; Ramella et al., 1989, 2001, 2002), or from projected (2-D) data (de Vaucouleurs, 1975a, b; Turner & Gott, 1976; Materne, 1978; de Filippis et al., 2000). All these catalogues are derived from different data sets and with different algorithms and are therefore affected by different biases favouring the detection of structures in a given richness range; biases induced by the topology of the data, by the limited size of the survey, by ambiguities in the selection criteria, etc. Shectman (1985) pioneered the field of automated cluster finding in optical surveys using peak-finding methods, which has been refined and modified in many later projects (Maddox et al., 1990; Dalton et al., 1992; Lumsden et al., 1992; Nichol et al., 2001a; Gal et al., 2002). Based on a model dependent approach, Postman et al. (1996) developed the matched filter technique, which has been widely used, with several variants (Kawasaki et al., 1998; Schuecker & Bohringer, 1998; Lobo et al., 2000), including the adaptive matched filter (Kepner et al., 1999). In addition, the availability of multiband high accuracy CCD data, allowed the implementation of several cluster-finding methods based on the use of galaxy colours (Gladders & Yee, 2000; Goto et al., 2002; Nichol et al., 2001b; Andreon, 2003). An independent approach relied on the Voronoi tessellation technique as a peak finder (Ramella et al., 2001; Kim et al., 2000) and a modified version, taking into account colours was implemented by Kim et al. (2002). More recently, other, more advanced pattern recognition tools such as Bayesian clustering (Murtagh et al., 2002), maximum likelihood (Cocco & Scaramella, 1999), and neural networks (Frattale Mascioli, Priv. Comm.) have been introduced.

Much less work has been done to detect poorer structures such as loose groups; two principal methods (and their successive elaborations) have been adopted. Turner & Gott (1976) presented the first tentative objective identification of groups as enhancements above a reliable threshold in the projected galaxy distribution. The ”Friends Of Friends” algorithm of Huchra & Geller (1982) generates a measure of correlation among galaxies and their neighbours, basing on their separation in the full 3-D space. A noticeable exception to the lack of low-richness catalogs has been the detection of compact groups, where several teams (de Carvalho & Djorgovski, 1995; Iovino et al., 1999, 2002) have proposed different approaches to their detection. For the determination of the MF, it is important to note that its derivation from the above cited catalogues is hindered by the fact that all of the above algorithms are optimised for the detection of either groups or clusters, and no systematic work has been done in matching their outcomes in the transition region between structures of low and high richness.

Here, we attempt the derivation of an accurate MF, starting from the galaxy catalogues extracted from DPOSS material.

The paper is structured as follows. In Sect. 2 we shortly summarise the properties of the Digitized Palomar Sky Survey (DPOSS) data (Djorgovski et al., 1998, 1999; Reid et al., 1991) used to derive the multiplicity function described in Sect. 5. In Sect. 3 we describe the algorithms used to detect groups (Sect. 3.1) and clusters (Sect. 3.2), while in Sect. 4 we discuss the simulations performed in order to evaluate the accuracy of the method, expressed in terms of completeness and fraction of spurious detections, and to evaluate the possible existence of systematic errors in the ranges of overlapping richness for the group and cluster finding procedures. Finally, in Sect. 6, we draw our conclusions. Through this paper we assume $`H_0=100`$ $`km`$ $`s^1Mpc^1`$.

## 2 The data

The data used in this paper were extracted from the DPOSS photographic plates (Djorgovski et al., 1998, 1999; Reid et al., 1991) using the SKICAT package (Weir et al., 1995a) which provides photometric, morphological and astrometric data for each detected object. SKICAT also provides a classification (Star/Galaxy) based on a classification tree (Weir et al., 1995b).
In DPOSS, the three photometric bands ($`J,F`$ and $`N`$) are individually calibrated to the Gunn system (Thuan & Gunn, 1976; Wade et al., 1979) by means of accurate CCD photometry of objects of intermediate luminosity, (to take into account the nonlinear response of the plates), with preferential targetting of galaxies. From the DPOSS data covering the selected regions, we extract, for each individual object: RA, Dec, total magnitude which best approximates the asymptotic magnitudes and the object classification.

DPOSS individual plate catalogues must be cleaned of spurious objects and artifacts (such as multiple detections coming from extended patchy objects, halos of bright stars, satellite tracks, etc.). In order to do so, we mask plate regions occupied by bright, extended and saturated objects which locally make object detection extremely unreliable. Subsequently, we matched catalogues obtained in each of the three photometric bands, by using the plate astrometric solution and by matching each object in one filter with the nearest objects in the two other filters (with a tolerance box of 7 arcsec, see Paolillo et al., 2001). Due to the different S/N ratios in the three bands, many objects had discordant star/galaxy classifications in catalogues obtained in the different bands. The number of such objects obviously increases at faint magnitudes (it needs to be stressed, however, that this problem is greatly reduced when a new training set for the classification is adopted, see Odewahn et al. (2002) for details). In order to exclude from our final catalogues the smallest number of true galaxies, we discard only the objects classified as stars in all three filters. Final catalogues were thus obtained for 10 DPOSS plates (see Tab. 1) covering a total area of $`300`$ sq. deg. spread at high galactic latitude ($`b>30\mathrm{deg}`$) (see Fig. 1), in order to reduce cosmic variance. Details on the photometric calibration of these particular plates can be found in Paolillo et al. (2001, 2003). We note that these calibrations are not the same as the general DPOSS calibrations described in Gal et al. (2002). Our catalogue of galaxies is limited in magnitude down to the Gunn $`r=20.5mag`$.

## 3 Detection of galaxy overdensities (groups and clusters)

Making an arbitrary choice, we use the term ”groups” to denote those galaxy aggregates which consist of less than 20 objects, and ”clusters” for all richer structures. This definition is comparable to that of Abell (1958), but in our case we set an implicit threshold on the magnitude difference between the brightest and the faintest objects in the same structure, given by the limiting magnitude.

### 3.1 The procedure for groups

In order to detect galaxy associations of low richness ($`N_{obj}<20`$), we have implemented a modified version of van Albada’s algorithm originally developed for binary systems (see Oosterloo, 1989; Soares et al., 1995).

Taking into account only the position and the apparent magnitude for each galaxy in our catalog, we first search for the nearest neighbour in a given magnitude range, and then estimate the probability that the two objects are physically related.

For the fore/background galaxies, the projected distribution is assumed to be Poissonian and the probability that the angular separation between a given galaxy and its nearest neighbour falls in the range $`\theta `$ and $`\theta +d\theta `$, is:

$$P_1(\theta )d\theta =\mathrm{exp}\left[\pi \theta ^2\rho \right]2\pi \theta \rho d\theta $$
(1)
where $`\rho `$ is the surface density of background galaxies in the immediate neighbourhood. In order to combine the angular separations of different pairs into a single distribution, the quantity $`x`$ is defined as the ratio between the observed value of the distance ($`\theta _1`$) to the nearest neighbour and the expected theoretical mean value $`<\theta _1>`$ given by Eq. 1:

$$\theta _1x\theta _1=\frac{x}{2\sqrt{\rho }}$$
(2)
The resulting frequency distribution of $`x`$:

$$p_1(x)dx=\mathrm{exp}\left[\frac{\pi x^2}{4}\right]\frac{\pi }{2}xdx$$
(3)
is then independent of the background density $`\rho `$.

The shape of the observed distribution, $`p_0(x)`$, and the Poisson distribution $`p_1(x)`$, for large $`x`$, are expected to be similar. If an excess is found in the observed distribution relative to the Poissonian expectation for small $`x`$ (see Fig. 2, lower panel), it is likely due to physical companions, which will tend to cluster at smaller distances than random projections.

Normalising the observed distribution to the Poisson distribution, we can use the excess $`p_0(x)p_1(x)`$, observed at small $`x`$, to define the probability $`p`$ that two galaxies, located at a certain $`x`$, are physically associated:

$$p1\frac{p_1(x)}{p_0(x)}$$
(4)
In this formalism, all galaxy pairs having $`p`$ higher than a given threshold value are considered to be physical companions.

Iteration of the above procedure allows us to estimate the probability that other companions of higher order (up to N $``$ 20) are physically related to the first object by comparing the observed distributions of higher order to the expected Poissonian distributions (normalised to the local density) for the second, third, etc. nearest neighbours ($`p_2(x)`$, $`p_3(x)`$, etc.).

Groups are then identified by associating all galaxies having probability $`p`$ higher than a given threshold value. Groups sharing one or more companions are finally merged into one single system. The total number of objects defines our richness for the groups.

To compute the quantity $`x`$ for every pair of galaxies, it is necessary to have an accurate estimate of the local galaxy density background $`\rho `$. To derive $`\rho `$, for each galaxy and within each magnitude interval, one first determines the distance to the $`i`$-th nearest neighbour $`\theta _i`$. The relation between $`\theta _i`$ and $`\rho `$ is given by the probability that the distance to the $`i`$-th nearest neighbour lies between $`\theta `$ and $`\theta +d\theta `$:

$$P_i(\theta )d\theta =\mathrm{exp}\left[\pi \theta ^2\rho \right]\frac{\left(\pi \theta ^2\rho \right)^{i1}}{\left(i1\right)!}2\pi \theta \rho d\theta $$
(5)
The mean expected value of $`\theta _i`$ is:

$$\theta _i=\frac{\mathrm{\Gamma }(i+0.5)}{[(i1)!\sqrt{(}\pi \rho )]}$$
(6)
The higher the chosen value of $`i`$ (i.e. for large distances to the $`i`$-th neighbours), the lower the probability of being affected by possible physical companions, which would lead to an overestimate of the local background. Furthermore the width of the distribution of the ratio between $`\theta _i`$ and its mean value ($`<\theta _i>`$) decreases with increasing $`i`$. Thus if $`i`$ is large enough, it is possible to obtain an accurate estimate of $`\rho `$ from Eq. 6 by replacing the expected value $`<\theta _i>`$ with the observed one $`\theta _i`$. On the other hand, $`i`$ must not be too large, otherwise too much of the small scale clustering would be lost, and a large area of the plate will be affected by border effects (distant companions of galaxies located near the border of the plate will not follow a Poissonian statistics and will be preferentially located towards the center of the plate).
The choice of the value of $`i`$ is therefore a compromise that has to be made by taking into account all ofthe above factors.

### 3.2 The procedure of cluster identification

Candidate clusters were identified following a procedure similar to that of Shectman (1985). The catalogue of galaxies is binned into equal-area square bins in the sky, generating a two dimensional map (density map) of the number density of galaxies. The bin size ($`1.2\mathrm{}`$) is chosen such that the mean number of galaxies per bin is $`1`$, in order not to degrade the spatial resolution. The resulting map (Fig. 3) exhibits irregular structures corresponding to the presence of overdensities emerging above the intrinsic fluctuations of the background distribution.

The large fluctuations existing in the distribution of background galaxies are due to the non-uniform background galaxy distribution. Once the density map has been created, the analysis of these maps poses similar problems to those of classical photometry, so we use S-Extractor (Bertin & Arnouts, 1996) for the detection of areas showing enhanced signal. S-Extractor is run on the density map searching for objects with a minimum detection area of 4 pixels above a global threshold of $`0.4`$ times the Poissonian background noise estimated from each plate using a background map. The evaluation of such background is a crucial step, strongly affecting the final richness estimate. The use of S-Extractor poses several problems (which cannot be trivially solved) since it is optimised to work on images with Gaussian statistics, while in density maps there are too few objects per bin, and they are distributed according to Poissonian statistics, thus making the background determination provided by S-Extractor unreliable. To circumvent this problem, we were forced to derive the background map in an alternative way. We first divide the original density map into sub-images of $`1^{}\times 1^{}`$, and then compute the Poissonian mean in each box, subsequently performing a fit with a 2-dimensional polynomial function of first order. We found a mean background density of 1640 per sq.deg. with a $`\sigma `$ of 148 galaxies. In this way, we remove those spatial frequencies higher (the clusters) than the mesh scale length. At the estimated typical redshift in our sample ($`z=0.10.2`$) this scale corresponds to a linear dimension of 9-15 Mpc. This map was then subtracted from the global frame before running the detection procedure.

The resulting density map was then smoothed in the detection step using S-Extractor with a Gaussian 2-D filter in order to match the cluster density profiles and, since we are searching structures with almost a Gaussian core, the filter width was chosen depending on the expected average apparent size for the cores ($``$ 250 kpc) of clusters in the redshift range ($`z=0.10.2`$) probed by our data. We stress that the choice of the otimal parameters strongly depends on the characteristics of the specific data sets and needs to be tuned on the simulations reproducing the behaviour of true catalogues.

The extracted parameters characterizing the detected overdensities are the density centroid in absolute equatorial coordinates (J2000), the isophotal area above the threshold, the S/N ratio of detection, and the number of objects inside the isophotal area, which we use to derive (after the background correction) our richness parameter for the clusters.

## 4 Outlines of the simulation

In order to test the limits of our group and cluster detection procedures, we performed simulations over a region having the same area and the galaxy counts as one POSS-II plate. In this way we could estimate the shortcomings of our procedure, such as the percentage of spurious detections and the percentage of lost objects; at the same time this helped in the fine tuning of the parameters of the detecting algorithms.

### 4.1 Simulation of the galaxy background

First we simulated the galaxy background assuming a uniform galaxy distribution. The number of simulated background galaxies is the average number of galaxies present in the DPOSS plates (approx $`50000`$ after excluding all the galaxies fainter than the limiting magnitude). To each background galaxy, a sky position, randomly extracted within the plate limits, and an apparent magnitude, distributed according to the observed galaxy counts, were assigned.

### 4.2 Simulation of galaxy groups

The number of groups to be simulated was extracted from the multiplicity function of Turner & Gott (1976).
We began by placing the principal galaxy of each group at random positions inside each field. Then, to each principal galaxy we assign an absolute magnitude and a redshift. Absolute magnitudes were extracted from a Schechter function with $`\mathrm{M}^{}=19.80`$ and $`\alpha =1.25`$ (Ramella et al., 1999), while the redshifts were assigned from the galaxy distribution observed in the Las Campanas redshift survey (Shectman et al., 1996). To each principal galaxy we then associate a number of secondary galaxies matching the multiplicity function mentioned above, each of these galaxies having the redshift of the corresponding principal galaxy. Taking into account the estimates provided in the literature, each simulated group was given a maximum standard dimension depending on its richness: a maximum radius of $`0.26\mathrm{Mpc}`$ for groups with $`N_{obj}18`$ members, while a maximum radius of $`0.55\mathrm{Mpc}`$ is used for groups with $`N_{obj}>18`$ members. All the secondary galaxies belonging to a group were then distributed inside the group volume, and each assigned an absolute magnitude generated from the same Schechter function as the brightest galaxies in the group. Finally, absolute magnitudes were re-transformed to apparent magnitudes by taking into account the cluster distance and the average $`k`$-corrections from Fukugita et al. (1995). The detection algorithm was then applied to the simulated plates in order to fine tune the algorithm parameters (threshold value of the probability $`p`$ and choice of the $`i\mathrm{th}`$ nearest neighbour to compute the background galaxy density). The results of the simulations may be summarised as follows: the group detection algorithm loses 28% of the simulated groups and produces 43% spurious detections.
Fig. 4 shows that, in spite of the high contamination level, the MF shape is statistically preserved: the simulated MF (filled circles) and the detected MF (empty squares) differ on average by a vertical offset, which we take into account to correct the final group MF.

In Fig. 5 (left panel), we show, as an example, the outcome of one typical simulation. The centers of the simulated (dark dots) and detected (empty triangles) groups are plotted; a circle is drawn when the two match.

### 4.3 Simulation of clusters

Cluster simulations were performed with same assumptions used for the groups, with some crucial differences. The number of simulated clusters of a given richness (ranging from 2 to 200 galaxies) in an area of 37.59 squared degrees (approximately the area of one DPOSS plate) was determined from a preliminary analysis performed on $`10`$ DPOSS plates. In a second step, a power law multiplicity function was used, with the slope taken from the preliminary multiplicity function. In this way we tried to take into account the total number of low richness objects, which could not be measured from our preliminary analysis.

The absolute magnitudes of the principal galaxies were extracted from a Gaussian distribution centered on $`22.99\pm 1.0\mathrm{mag}`$ (Schneider et al., 1983), while those for the secondary galaxies were extracted from the luminosity function of Paolillo et al. (2001). To take into account the richness dependence of the cluster dimensions, we arbitrarily adopted a core radius ($`\sigma `$ of the Gaussian profile) of $`0.5\mathrm{Mpc}`$ for clusters with $`30`$ members, while a core radius of $`1.0\mathrm{Mpc}`$ was used for clusters with $`>30`$ members. Although these values may appear somewhat high, the adoption of smaller values for the core radius would only make the detection easier and therefore the whole procedure more reliable. As with the groups, the detection algorithm was applied to a large number of simulated plates to test the algorithm performance as a function of the properties of the objects to be detected.

In Fig. 6, we plot, for a typical simulated plate, the assigned richness vs. the assigned core radius of the simulated clusters (open circles) and mark with a cross the clusters retrieved by the algorithm. Clusters with a very shallow profile or which are poor are preferentially lost.

The dependence of the algorithm efficiency on the richness is shown in Fig. 7, where we plot the number of simulated (continuous line) and retrieved (dash shaded area) clusters in the typical plate area. All but two of the clusters having $`N_{obj}>35`$ are retrieved. In the range of richness $`25<N_{obj}<35`$, 80% of the clusters are retrieved. Considering that a cluster belonging to the Abell richness class 0 ($`3049`$ members in a range of two magnitudes) has $`N_{obj}>30`$ ($`N_{obj}`$ includes the cluster galaxies in a range of at least four magnitudes), we are complete up to $`z=0.2`$ at least for all the Abell richness classes. Fig. 7 also shows that spurious detections (dot shaded area) are absent in the richness range where the algorithm works with the highest efficiency, and occur only in the range where the group finder is to be used.

As already mentioned, the estimate of clusters richness is given by the number of objects within the detection isophote (isodensity counts). We wish to stress that this definition of richness depends on the redshift of the detected structure. The quality of the richness estimate has been tested using our simulations. In Fig. 8, the points follows bisector of the diagram (a bit shifted towards the upper half part of the plot), with a scatter in richness of $`10`$ galaxies (which is consistent with the background fluctuations). The small shift indicates an understimation of the retrieved richnesses. We are comparing the number of the galaxies put in a synthetic circular aperture (the simulated) with the richness in the isodensity irregular countours, as it is measured in the real case: in this way some galaxies are missed. If we use circular apertures of the cluster size (which are known in the simulations but not in the actual observations), the shift disappears. Points in the lower right part of the plot are due to overlapping clusters, for which (in the absence of a deblending procedure) the richness will obviously be overestimated.

## 5 The conjoined groups/clusters multiplicity function

Fig. 9 summarizes our main results. We plot the MF, defined as the number of groups or clusters per unit area and per unit of estimated richness (the groups/clusters richness is defined respectively in Sect. 3.1 and Sect. 3.2). For the clusters, the bin grows exponentially as $`2^{n/2}`$, in order to keep the $`S/N`$ ratio almost constant along the richness axis. For the groups, the bin was instead set equal to 1. In order to exclude the structures detected in the redshift range where our magnitude limited catalogue is incomplete, only clusters and groups where the brightest galaxy has $`m<16.5`$ (in Gunn $`r`$) were selected. Assuming that brightest galaxies may be used as standard candles, our selection in magnitude implies that $`z<0.2`$.

The procedures described above were applied separately for groups and clusters, obtaining two different multiplicity functions (marked with different symbols in Fig. 9). These MFs appear to define a common relation, without the need for any offsets or normalisations. We emphasize that a minor correction for completeness was applied only to the last point of the clusters MF. To correct the group’s MF for contamination by spurious detections (see Sect. 4.2), a global shift derived from the simulations was also applied. Only the Poissonian statistical fluctuations have been taken into account in the error estimate. For the high richness clusters, the error on the richness estimate is negligible with respect to the bin width. The error becomes relevant only in the same richness range where incompleteness is also significant.

In Fig. 10 we compare our results with a MF extracted by us from the USGC catalog of groups (Ramella et al., 2002). We adopt the same representation scheme for the two data set. Normalisation to the same volume was applied to the USGC groups, assuming a uniform distribution of objects in redshift both for our sample and the USGC sample. It is important to note that the two catalogs were derived in totally different ways. The USGC is generated using spectroscopic redshifts by a percolation method implemented by Ramella et al. (1997) for group detection, which is designed to reduce the risk of false detections introduced by chance projections.

The agreement between these two MFs (see Fig. 10), derived under totally different assumptions and using independent data sets, is due to similar biases affecting the estimated richnesses for both samples. For low $`N_{obj}`$ structures (groups) the similarity is apparent; in both cases, the methods count individual objects fulfilling the respective membership criteria but with secondary members having magnitudes falling within similar (i.e. $``$ four magnitudes) ranges with respect to the primary galaxy. For the clusters, instead, the different depths sampled by the two data sets, when compared to the different limiting magnitudes of the samples themselves, indicates that both methods sample very similar intervals of the cluster’s luminosity function.

## 6 Summary and discussion

We have implemented two algorithms for the detection of galaxy associations, one for groups and one for clusters. The former is a modified version of van Albada’s procedure to detect galaxy pairs, while the latter the Shectman (1985) approach, which uses a peak-finding procedure on a density map obtained from the galaxy catalogue.

We evaluated the performance of these methods via extensive simulations, which show that the group algorithm is reliable up to richness $`20`$, and the cluster algorithm is reliable at richnesses above $`20`$ galaxies. The two algorithms were then applied to a $`300`$ square degree field extracted from DPOSS data (see Sect. 2). The resulting MFs show a remarkable internal consistency from the two procedures which produce independent MFs for groups and clusters, matching with no need for normalisation. Additionally, the MF derived using our technique on the 3-D based catalogues of Ramella et al. (2002) agrees with the MF derived from the projected DPOSS data. The final combined MF is well fit by a power-law of slope $`\alpha =2.08\pm 0.07`$. The correlation coefficient on the log-log scale is $`0.98`$.

The data set we used to determine the MF samples a volume \[$`300deg^2`$, $`z<0.2`$\] which is slightly smaller than that explored by Bachall et al. (2002) \[$`400deg^2`$, $`z=0.10.2`$\]. The total number of detected structures for $`N>10`$ in the Bachall et al. (2002) and in our sample is respectively $`300`$ and $`370`$. In a forthcoming paper we shall analyze the cosmological implications of the derived MF.

###### Acknowledgements.

The authors wish to thank Marisa Girardi and Michail Sazhin for useful comments and stimulating discussions.