text:
# A Toy Model for the Magnetic Connection between a Black Hole and a Disk

## I Introduction

A magnetic field connecting a Kerr black hole to a disk rotating around it can transfer energy and angular momentum between the black hole and the disk (bla99 ; li00a ; li00b ; li00c  and references therein). When the black hole rotates faster than the disk, energy and angular momentum are extracted from the black hole and transfered to the disk. If the magnetic field is weak, the energy transfered to the disk by the magnetic connection can be dissipated and radiated away by the disk li00a ; li00b ; li00c . Thus, in the presence of magnetic connection, the rotational energy of the black hole provides an energy source for disk radiation in addition to disk accretion nov73 ; pag74 ; tho74 . In fact, a disk magnetically coupled to a rapidly rotating Kerr black hole can radiate without the need of accretion li00a ; li00b ; li00c . In such a case, all the power of the disk comes from the rotational energy of the black hole.

Due to the lack of knowledge about the topology and distribution of the magnetic field connecting the black hole to the disk, a detailed model of the magnetic connection between a black hole and a disk does not exist yet, except the simplest and ad hoc model where the magnetic field is assumed to connect the black hole to the disk close to the inner boundary li00a ; li00b , or at a single radius beyond the inner boundary li00c . In this paper, we consider a toy model for the magnetic connection: A single toroidal electric current flows around a Kerr black hole in the equatorial plane which is coincident with the plane of a thin Keplerian disk rotating around the black hole. The toroidal current generates a poloidal magnetic field which connects the black hole to the disk. Due to the relative rotation between the black hole and the disk and the fact that the horizon of a Kerr black hole is a conductor with a finite electric resistivity zna78 ; dam78 ; car79  and a plasma disk is a perfect conductor, an electromotive forces (EMF) is induced which in turn generates a poloidal electric current flowing through the black hole and the disk. The poloidal current produces a torque and a power on the disk, which makes the disk shine without the need of accretion. We will calculate the power and the resulting radiation flux of the disk with the assumption that there is no accretion and the configuration is steady and axisymmetric, and compare the results with that of a standard accretion disk.

A similar model has been used by Li li00d  to test the efficiency of the Blandford-Znajek mechanism bla77 ; mac82 ; phi83 . By assuming a toroidal current continuously distributed in a thin Keplerian disk, Li li00d  has shown that the power of the disk always dominates over the power of the black hole, which confirms the early speculation of Blandford and Znajek bla77  and the recent arguments of Livio, Ogilvie, and Pringle liv99 . In this paper we will find that, an electric current flowing around a black hole in the equatorial plane naturally generates a magnetic field connecting the black hole to its disk. This implies that the existence of magnetic connection between a black hole and a disk is a very natural assumption in view of physics. Essentially the magnetic connection between a black hole and a disk is a variant of the Blandford-Znajek mechanism. However, it is simpler and better defined than the Blandford-Znajek mechanism since a unknown remote load is not required li00c . Though the model presented in this paper is so simple that it cannot be directly applied to real astronomy, it helps us in understanding some physical aspects of the magnetic connection.

Throughout the paper we use the geometrized units $`G=c=1`$ and the Boyer-Lindquist coordinates $`(t,r,\theta ,\varphi )`$ mis73 ; wal84 .

## II The magnetic flux and the magnetic connection

Assuming a neutral toroidal electric current $`I`$ flows on a circle of $`r=r^{}`$ in the equatorial plane of a Kerr black hole, then the four-current vector is

$`J^a={\displaystyle \frac{I}{r}}\left({\displaystyle \frac{\mathrm{\Delta }}{A}}\right)^{1/2}\left({\displaystyle \frac{}{\varphi }}\right)^a\delta \left(rr^{}\right)\delta (\mathrm{cos}\theta ),`$ (1)
where $`\mathrm{\Delta }r^22Mr+a^2`$, $`A(r^2+a^2)^2\mathrm{\Delta }a^2\mathrm{sin}^2\theta `$, $`M`$ is the mass of the black hole, $`Ma`$ is the angular momentum of the black hole ($`a^2M^2`$), and $`\delta (x)`$ is the Dirac $`\delta `$-function. The magnetic flux through a surface bounded by a circle with $`r=\mathrm{constant}`$ and $`\theta =\mathrm{constant}`$, produced by such an electric current, is

$`\mathrm{\Psi }(r,\theta ;r^{})=2\pi A_\varphi (r,\theta ;r^{}),`$ (2)
where $`A_\varphi `$ , the toroidal component of the electric vector potential, is given by zna78a ; lin79

$`A_\varphi `$ $`=`$ $`2{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}\{\alpha _l^r[ra\mathrm{sin}^2\theta {\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Sigma }}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}P_l^{}(u)P_l(\mathrm{cos}\theta )a\mathrm{sin}^2\theta \mathrm{cos}\theta {\displaystyle \frac{r^2+a^2}{\mathrm{\Sigma }}}P_l(u)P_l^{}(\mathrm{cos}\theta )]`$ (3)
$`+\alpha _l^i[a^2\mathrm{sin}^2\theta \mathrm{cos}\theta {\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Sigma }}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}P_l^{}(u)P_l(\mathrm{cos}\theta )r\mathrm{sin}^2\theta {\displaystyle \frac{r^2+a^2}{\mathrm{\Sigma }}}P_l(u)P_l^{}(\mathrm{cos}\theta )`$
$`+{\displaystyle \frac{\mathrm{\Delta }\mathrm{sin}^2\theta }{l(l+1)}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}P_l^{}(u)P_l^{}(\mathrm{cos}\theta )]\}`$
$`+2{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}\{\beta _l^r[ra\mathrm{sin}^2\theta {\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Sigma }}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}Q_l^{}(u)P_l(\mathrm{cos}\theta )a\mathrm{sin}^2\theta \mathrm{cos}\theta {\displaystyle \frac{r^2+a^2}{\mathrm{\Sigma }}}Q_l(u)P_l^{}(\mathrm{cos}\theta )]`$
$`+\beta _l^i[a^2\mathrm{sin}^2\theta \mathrm{cos}\theta {\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Sigma }}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}Q_l^{}(u)P_l(\mathrm{cos}\theta )r\mathrm{sin}^2\theta {\displaystyle \frac{r^2+a^2}{\mathrm{\Sigma }}}Q_l(u)P_l^{}(\mathrm{cos}\theta )`$
$`+{\displaystyle \frac{\mathrm{\Delta }\mathrm{sin}^2\theta }{l(l+1)}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}Q_l^{}(u)P_l^{}(\mathrm{cos}\theta )]\},`$
where $`u(rM)/\sqrt{M^2a^2}`$, $`P_l(z)`$ and $`Q_l(z)`$ are Legendre functions, $`P_l^{}(z)dP_l(z)/dz`$, $`Q_l^{}(z)dQ_l(z)/dz`$, and the coefficients $`\alpha _l^r`$, $`\alpha _l^i`$, $`\beta _l^r`$, and $`\beta _l^i`$ are respectively
(1) for $`r<r^{}`$, $`\beta _l^r=\beta _l^i=0`$ for all $`l`$; but

$`\alpha _l^r{\displaystyle \frac{\pi (2l+1)I}{l(l+1)(M^2a^2)}}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{A^{}}}\right)^{1/2}\mathrm{\Delta }^{}aP_l(0)Q_l^{}(u^{}),`$ (4)
$`\alpha _l^i`$ $``$ $`{\displaystyle \frac{\pi (2l+1)I}{l(l+1)\sqrt{M^2a^2}}}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{A^{}}}\right)^{1/2}[(r_{}^{}{}_{}{}^{2}+a^2)P_l^{}(0)`$ (5)
$`\times Q_l(u^{})+{\displaystyle \frac{r^{}\mathrm{\Delta }^{}}{l(l+1)}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}P_l^{}(0)Q_l^{}(u^{})];`$
(2) for $`r>r^{}`$, $`\alpha _l^r=\alpha _l^i=0`$ for all $`l`$; but

$`\beta _l^r{\displaystyle \frac{\pi (2l+1)I}{l(l+1)(M^2a^2)}}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{A^{}}}\right)^{1/2}\mathrm{\Delta }^{}aP_l(0)P_l^{}(u^{}),`$ (6)
$`\beta _l^i`$ $``$ $`{\displaystyle \frac{\pi (2l+1)I}{l(l+1)\sqrt{M^2a^2}}}\left({\displaystyle \frac{\mathrm{\Delta }^{}}{A^{}}}\right)^{1/2}[(r_{}^{}{}_{}{}^{2}+a^2)P_l^{}(0)`$ (7)
$`\times P_l(u^{})+{\displaystyle \frac{r^{}\mathrm{\Delta }^{}}{l(l+1)}}{\displaystyle \frac{1}{\sqrt{M^2a^2}}}P_l^{}(0)P_l^{}(u^{})];`$
where $`\mathrm{\Delta }^{}=\mathrm{\Delta }(r=r^{})`$ and $`A^{}=A(r=r^{},\theta =\pi /2)`$.

The magnetic flux through the hemi-sphere of the horizon of the black hole is $`\mathrm{\Psi }(r_H,\pi /2;r^{})`$, where $`r_H=M+\sqrt{M^2a^2}`$ is the radius of the horizon. The magnetic flux through a surface bounded by a circle $`r=\mathrm{constant}`$ in the equatorial plane is $`\mathrm{\Psi }(r,\pi /2;r^{})`$. The magnetic field lines generated by the toroidal current are closed loops around the circle $`r=r^{}`$, except the field line along the axis $`\theta =0`$. Consider the field lines leaving the horizon at a polar angle $`\theta _1`$, they will intersect the equatorial plane at a circle of radius $`r_1`$ where $`r_1`$ is given by $`\mathrm{\Psi }(r_1,\pi /2;r^{})=\mathrm{\Psi }(r_H,\theta _1;r^{})`$. Clearly, if the inner edge of the disk lying in the equatorial plane is inside the circle $`r=r_1`$ and the outer edge of the disk is outside the circle $`r=r_1`$, a magnetic connection must exist between the black hole and the disk.

Let us consider a thin Keplerian disk in the equatorial plane of a Kerr black hole nov73 ; pag74 ; tho74 . The angular velocity of the disk is

$`\mathrm{\Omega }_D(r)=\left({\displaystyle \frac{M}{r^3}}\right)^{1/2}{\displaystyle \frac{1}{1+a\left(\frac{M}{r^3}\right)^{1/2}}}.`$ (8)
The outer edge of the disk is at $`r=\mathrm{}`$. The inner edge of the disk is at the marginally stable orbit

$`r_{ms}=M\left\{3+z_2\left[(3z_1)(3+z_1+2z_2)\right]^{1/2}\right\}`$ (9)
where

$`z_1`$ $``$ $`1+\left(1{\displaystyle \frac{a^2}{M^2}}\right)^{1/3}`$ (10a)
$`\times \left[\left(1+{\displaystyle \frac{a}{M}}\right)^{1/3}+\left(1{\displaystyle \frac{a}{M}}\right)^{1/3}\right],`$
$`z_2`$ $``$ $`\left(z_1^2+3{\displaystyle \frac{a^2}{M^2}}\right)^{1/2}.`$ (10b)
Due to the existence of a gap between the horizon of the black hole and the inner boundary of the the disk, not all magnetic field lines connect the black hole to the disk. If we define a flux surface to be a surface of constant magnetic flux, there is an innermost flux surface which connects the black hole to the disk. Assuming the toroidal electric current is located at the inner boundary of the disk: $`r^{}=r_{ms}`$, then the innermost flux surface touches the horizon of the black hole at $`\theta =\pi /2`$ and intersects the disk at a circle whose radius $`r_{in}`$ is determined by

$`\mathrm{\Psi }(r_{in},{\displaystyle \frac{\pi }{2}};r_{ms})=\mathrm{\Psi }(r_H,{\displaystyle \frac{\pi }{2}};r_{ms}),r_{in}>r_{ms}.`$ (11)
Since the magnetic field lines in the gap between the black hole and the disk must close themselves on the disk, $`r_{in}`$ must be greater than $`r_{ms}`$. Clearly, $`r_{in}/M`$ depends only on $`a/M`$ when $`r^{}=r_{ms}`$ is fixed. In Fig. 1 we plot $`r_{in}/M`$ as a function of $`a/M`$. For comparison, we also plot $`r_{ms}/M`$ and $`r_H/M`$ in Fig. 1. We see that, as $`a/M`$ increases, $`r_{in}`$ gets closer to $`r_{ms}`$. But, for $`a/M<1`$, we always have $`r_{in}>r_{ms}`$. The magnetic field lines between $`r_{ms}`$ and $`r_{in}`$ close themselves in the gap region between $`r_H`$ and $`r_{ms}`$, they do not touch the black hole horizon.

In order for the magnetic field to connect the black hole to the whole disk from $`r=r_{ms}`$ to $`r=\mathrm{}`$, the toroidal current must be inside the inner edge of the disk but outside the horizon of the black hole. Assume $`r^{}=r_0`$ when the inner most flux surface touches the disk at $`r=r_{ms}`$ and the horizon of the black hole at $`\theta =\pi /3`$ (and $`\theta =2\pi /3`$ due to the reflection symmetry about the equatorial plane). Then, $`r_0`$ is determined by

$`\mathrm{\Psi }(r_{ms},{\displaystyle \frac{\pi }{2}};r_0)=\mathrm{\Psi }(r_H,{\displaystyle \frac{\pi }{3}};r_0),r_H<r_0<r_{ms}.`$ (12)
Like $`r_{in}/M`$, $`r_0/M`$ also depends only on $`a/M`$. In Fig. 2 we plot $`r_0/M`$ as a function of $`a/M`$. For comparison, we also plot $`r_{ms}/M`$ and $`r_H/M`$ in Fig. 2.

Examples of the topological structure of the magnetic field produced by a toroidal electric current flowing in the equatorial plane of a Kerr black hole are shown in Fig. 3 and Fig. 4.

## III The power produced by the magnetic connection

In the presence of a magnetic field a Kerr black hole behaves like a conductor with a surface resistivity $`R_H=4\pi 377`$ Ohms zna78 ; dam78 ; car79 . A plasma disk can be treated as a perfect conductor with a vanishing resistivity. So, the rotation of the black hole and the disk induces an EMF on the horizon of the black hole and an EMF on the disk mac82 ; li00d . For a bunch of magnetic field lines connecting the black hole to the disk within radii $`r`$$`r+\mathrm{\Delta }r`$, where $`\mathrm{\Delta }rr`$, the induced EMFs on the horizon of the black hole and the disk are respectively

$`\mathrm{\Delta }_H={\displaystyle \frac{1}{2\pi }}\mathrm{\Omega }_H\mathrm{\Delta }\mathrm{\Psi },\mathrm{\Delta }_D={\displaystyle \frac{1}{2\pi }}\mathrm{\Omega }_D\mathrm{\Delta }\mathrm{\Psi },`$ (13)
where $`\mathrm{\Delta }\mathrm{\Psi }`$ is the magnetic flux connecting the black hole to the disk, $`\mathrm{\Omega }_H`$ is the angular velocity of the black hole

$`\mathrm{\Omega }_H={\displaystyle \frac{a}{2Mr_H}},`$ (14)
which is constant over the horizon of the black hole. For a thin Keplerian disk the angular velocity of the disk is given by Eq. (8).

The black hole and the disk form a closed electric circuit, an electric current flows along the magnetic field lines connecting the black hole to the disk and closes itself inside the disk and the black hole. We assume that there is a thin plasma corona around the black hole and the disk, which can carry the poloidal current flowing between the black hole and the disk bla77 . In the magnetosphere the resistivity along the magnetic field lines is negligible, while the resistivity perpendicular to the magnetic field lines is large. Thus the electric current flows along the magnetic field lines in the corona without dissipation. Suppose the disk and the black hole rotates in the same direction, then the EMF of the black hole, $`\mathrm{\Delta }_H`$, and the EMF of the disk, $`\mathrm{\Delta }_D`$, have opposite signs. The direction of the electric current, and thus the direction of the transfer of energy and angular momentum, is determined by the sign of $`\mathrm{\Delta }_H+\mathrm{\Delta }_D`$. If $`\mathrm{\Delta }_H+\mathrm{\Delta }_D>0`$, the EMF of the black hole dominates the EMF of the disk, so the black hole “charges” the disk: energy and angular momentum are transfered from the black hole to the disk. On the other hand, if $`\mathrm{\Delta }_H+\mathrm{\Delta }_D<0`$, the EMF of the disk dominates the EMF of the black hole, so the disk “charges” the black hole, energy and angular momentum are transfered from the disk to the black hole li00b .

An infinite number of adjacent infinitesimal poloidal electric current loops flow between the black hole and the disk along the magnetic field lines connecting them. Each infinitesimal current loop produces an infinitesimal torque and an infinitesimal power on the disk li00a

$`\mathrm{\Delta }T_{HD}=\left({\displaystyle \frac{\mathrm{\Delta }\mathrm{\Psi }}{2\pi }}\right)^2{\displaystyle \frac{\mathrm{\Omega }_H\mathrm{\Omega }_D}{\mathrm{\Delta }Z_H}},`$ (15)
$`\mathrm{\Delta }P_{HD}=\mathrm{\Delta }T_{HD}\mathrm{\Omega }_D,`$ (16)
where $`\mathrm{\Delta }Z_H`$ is the resistance of the black hole associated with the infinitesimal current. In Eq. (16) we have used the fact that the disk is perfectly conducting, so the magnetic field lines are frozen in and corotate with the disk. The summation of all $`\mathrm{\Delta }T_{HD}`$ gives the the total torque on the disk, the summation of all $`\mathrm{\Delta }P_{HD}`$ gives the the total power on the disk. In the limit $`\mathrm{\Delta }r0`$ and $`\mathrm{\Delta }\mathrm{\Psi }0`$, the summation is replaced by integration. Then, after taking into account that a disk has two surfaces, we obtain the total torque and the total power on the disk

$`T_{HD}=4\pi {\displaystyle _{r_a}^{\mathrm{}}}Hr𝑑r,P_{HD}=4\pi {\displaystyle _{r_a}^{\mathrm{}}}\mathrm{\Omega }Hr𝑑r,`$ (17)
where $`r_a`$ is the radius in the disk where the magnetic connection starts ($`r_a=r_{in}`$ for $`r^{}=r_{ms}`$; $`r_a=r_{ms}`$ for $`r^{}=r_0`$),

$`H{\displaystyle \frac{1}{8\pi ^2r}}\left({\displaystyle \frac{d\mathrm{\Psi }}{dr}}\right)^2{\displaystyle \frac{\mathrm{\Omega }_H\mathrm{\Omega }_D}{dZ_H/dr}},`$ (18)
where we have treated the resistance of the black hole, $`Z_H`$, as a function of the radius of the disk, which is defined by a map from the horizon of the black hole to the surface of the disk given by the magnetic field lines connecting them.

For a Kerr black hole we have mac82

$`{\displaystyle \frac{dZ_H}{d\theta }}={\displaystyle \frac{R_H}{2\pi }}{\displaystyle \frac{r_H^2+a^2\mathrm{cos}^2\theta }{\left(r_H^2+a^2\right)\mathrm{sin}\theta }},`$ (19)
where $`\theta `$ is the polar angle coordinate on the horizon. Then, $`dZ_H/dr`$ can be calculated through

$`{\displaystyle \frac{dZ_H}{dr}}={\displaystyle \frac{dZ_H}{d\theta }}{\displaystyle \frac{d\theta }{dr}},`$ (20)
where $`\theta =\theta (r)`$ is a map between the $`\theta `$ coordinate on the horizon and the $`r`$ coordinate on the disk, which is induced by the magnetic field lines connecting the disk to the horizon

$`\mathrm{\Psi }(r_H,\theta ;r^{})=\mathrm{\Psi }(r,{\displaystyle \frac{\pi }{2}};r^{}).`$ (21)
Since $`d\theta /dr<0`$ and $`dZ_H/d\theta >0`$, we have $`dZ_H/dr<0`$.

For a given toroidal current loop at $`r=r^{}`$, we can calculate the magnetic flux with Eqs. (27), determine the innermost flux surface and thus $`r_a`$, determine the map defined by Eq. (21). Then, with Eq. (17) and Eq. (18), we can calculate the power of energy transfered by the magnetic connection. Let us define a characteristic magnetic field strength by

$`B_0{\displaystyle \frac{2I}{M}},`$ (22)
then we have $`\mathrm{\Psi }=\overline{\mathrm{\Psi }}B_0M^2`$, $`T_{HD}=\overline{T}B_0^2M^3`$, $`P_{HD}=\overline{P}B_0^2M^2`$, where the dimensionless quantities $`\overline{\mathrm{\Psi }}`$, $`\overline{T}`$, and $`\overline{P}`$ depend only on $`a/M`$ and $`r^{}/M`$. We have calculated $`\overline{P}`$ for two cases: (1) $`r^{}=r_{ms}`$, i.e., the toroidal current flows on the inner boundary of the disk; the radius where the innermost flux surface touches the disk ($`r_{in}`$) is shown in Fig. 1; (2) $`r^{}=r_0`$ where $`r_0`$ is defined by Eq. (12), so that the magnetic field connects the portions of the horizon of the black hole with $`0<\theta <\pi /3`$ and $`2\pi /3<\theta <\pi `$ to the whole disk from $`r=r_{ms}`$ to $`r=\mathrm{}`$; the values of $`r_0/M`$ as a function of $`a/M`$ are plotted in Fig. 2. The results for $`\overline{P}=P/P_0`$, where $`P_0B_0^2M^2`$, are shown in Fig. 5 and Fig. 6, respectively for the two cases.

In Fig. 5, the toroidal electric current flows on the marginally stable orbit (i.e. the inner edge of the disk). The poloidal magnetic field generated by such a current connects the horizon of the black hole to the regions of the disk with $`r>r_{in}`$, where $`r_{in}>r_{ms}`$ is defined by Eq. (11). Steady state solutions exist only for $`a/M0.01544`$. The power produced by the magnetic connection peaks at $`a/M=0.986`$ and drops to zero as $`a/M`$ approaches $`1`$, which is a magnification of the fact that an extreme Kerr black hole expels magnetic fields bic85 . If the black hole is spun down from $`a/M=0.998`$ to $`a/M=0.01544`$, about $`0.03Mc^2`$ energy is extracted from the black hole. In Fig. 6, the toroidal electric current flows on a circle of radius $`r_0`$ defined by Eq. (12). The poloidal magnetic field generated by such a current connects the portions of the horizon of the black hole with $`0<\theta <\pi /3`$ and $`2\pi /3<\theta <\pi `$ to the whole disk from $`r=r_{ms}`$ to $`r=\mathrm{}`$. Steady solutions exist only for $`a/M0.3594`$. The power produced by the magnetic connection peaks at $`a/M=0.961`$ and drops to zero as $`a/M`$ approaches $`1`$. If the black hole is spun down from $`a/M=0.998`$ to $`a/M=0.3594`$, about $`0.09Mc^2`$ energy is extracted from the black hole.

## IV The radiation flux of the disk

The energy pumped into the disk by the magnetic connection can be dissipated and radiated away by the disk. In a steady state without accretion, the radiation flux of the disk, which is the energy radiated per unit time and per unit area by the disk as measured by an observer corotating with the disk, is li00c

$`F=\{\begin{array}{cc}\frac{1}{r}\left(\frac{d\mathrm{\Omega }_D}{dr}\right)\left(E^+\mathrm{\Omega }_DL^+\right)^2_{r_a}^r\left(E^+\mathrm{\Omega }_DL^+\right)Hr𝑑r,\hfill & \text{if }r>r_a\hfill \\ 0,\hfill & \text{if }r<r_a\hfill \end{array},`$ (25)
where $`H`$ is defined by Eq. (18), $`E^+`$ and $`L^+`$ are the specific energy and the specific angular momentum of a particle in the disk, respectively. For a thin Keplerian disk, we have nov73 ; pag74

$`E^+={\displaystyle \frac{12\frac{M}{r}+a\left(\frac{M}{r^3}\right)^{1/2}}{\sqrt{13\frac{M}{r}+2a\left(\frac{M}{r^3}\right)^{1/2}}}},`$ (26)
$`L^+=\left(Mr\right)^{1/2}{\displaystyle \frac{12a\left(\frac{M}{r^3}\right)^{1/2}+\frac{a^2}{r^2}}{\sqrt{13\frac{M}{r}+2a\left(\frac{M}{r^3}\right)^{1/2}}}},`$ (27)
and

$`E^+\mathrm{\Omega }_DL^+={\displaystyle \frac{\sqrt{13\frac{M}{r}+2a\left(\frac{M}{r^3}\right)^{1/2}}}{1+a\left(\frac{M}{r^3}\right)^{1/2}}}.`$ (28)
Since $`r_{ms}r_a`$, $`F`$ is always zero at the inner boundary of the disk. As $`r\mathrm{}`$, $`H`$ declines quickly and the integration in Eq. (25) converges. Thus, as $`r\mathrm{}`$, $`F`$ approaches $`r^{3.5}`$.

As in Sec. III, we can calculate the magnetic flux with Eqs. (27) for a given toroidal current loop at $`r=r^{}`$, then determine $`r_a`$ and the map defined by Eq. (21). Then, with Eqs. (1820), Eq. (25), and Eq. (28), we can calculate the radiation flux of the disk. Using the characteristic magnetic field strength defined by Eq. (22), we have $`F=\overline{F}F_0`$ where $`F_0B_0^2`$, the dimensionless quantity $`\overline{F}`$ depends on $`a/M`$ and $`r^{}/M`$ only. We have calculated $`\overline{F}`$ for the following two cases: (1) $`r^{}=r_{ms}`$, $`a/M=0.99`$; (2) $`r^{}=r_0`$, $`a/M=0.99`$. The results are shown in Fig. 7 and Fig. 8, respectively. For comparison, we have also shown the radiation flux of a standard thin Keplerian disk (see Appendix A) around the same Kerr black hole in Fig. 7 and Fig. 8.

From Fig. 7 and Fig. 8, we see that the radiation flux produced by the magnetic connection is very different from that produced by the standard accretion: They peak at different radii, and at large radii the radiation flux produced by the magnetic connection declines faster than that produced by the standard accretion. The position of the peak of the radiation flux produced by the magnetic connection depends on the position of the toroidal electric current that generates the poloidal magnetic field. To quantify the variation speed of the radiation flux with radius, we have calculated the emissivity defined by $`\alpha d\mathrm{ln}F/d\mathrm{ln}r`$, the results are shown in Fig. 9. We see that, the magnetic connection can produce a very large emissivity index compared to the standard accretion, especially when the magnetic field gets more concentrated toward the center of the disk. As $`r\mathrm{}`$, the emissivity index of the magnetic connection approaches $`3.5`$, while the emissivity index of accretion approaches $`3`$.

## V Conclusions

In this paper we have constructed a toy model for the magnetic connection between a black hole and a disk. Assuming a toroidal electric current flows on a circle in the equatorial plane of a Kerr black hole, we have calculated the generated poloidal magnetic field which connects the horizon of the black hole to a thin Keplerian disk around the black hole, the power produced by the magnetic connection, and the radiation flux of the disk with the assumption of axisymmetry and steady state. We have also compared the results with that of a standard accretion disk.

We have particularly considered two cases: (1) The toroidal electric current flows on the marginally stable orbit where the inner boundary of the disk is. The magnetic field produced by such a toroidal current connects the horizon of the black hole to the disk from $`r=r_{in}`$ to $`r=\mathrm{}`$, where $`r_{in}>r_{ms}`$ is defined by Eq. (11). (2) The toroidal electric current flows on a circle of radius $`r_0`$ inside the marginally stable orbit, where $`r_0<r_{ms}`$ is defined by Eq. (12). The magnetic field produced by such a toroidal current connects the portions of the horizon of the black hole with $`0<\theta <\pi /3`$ and $`2\pi /3<\theta <\pi `$ to the whole disk from $`r=r_{ms}`$ to $`r=\mathrm{}`$.

We see that, a magnetic connection between a black hole and a disk is naturally produced by a toroidal current flowing around the black hole in the equatorial plane. The magnetic connection produces a power on the disk without need of accretion, the rotational energy of the black hole provides the energy source. The radiation flux of the disk produced by the magnetic connection is very different from that produced by accretion. Though the position of the peak of the radiation flux depends on the position of the toroidal electric current, at large radii the radiation flux of magnetic connection always declines faster than that of accretion. Over a large range of radii the emissivity index of magnetic connection is bigger than that of accretion. At large radii the emissivity index of magnetic connection approaches $`3.5`$ from above while the emissivity index of accretion approaches $`3`$ from below.

###### Acknowledgements.

I am very grateful to Bohdan Paczyński for stimulating discussions. This work was supported by a Harold W. Dodds Fellowship of Princeton University.

## Appendix A Radiation Flux of a Standard Accretion Disk

The radiation flux of a standard thin Keplerian disk around a Kerr black hole is nov73 ; pag74

$`F_{acc}={\displaystyle \frac{1}{4\pi r}}\dot{M}_Df,`$ (29)
where $`\dot{M}_D`$ is the mass accretion rate and

$`f`$ $`=`$ $`{\displaystyle \frac{3}{2M}}{\displaystyle \frac{1}{x^2\left(x^33x+2s\right)}}[xx_0{\displaystyle \frac{3}{2}}s\mathrm{ln}\left({\displaystyle \frac{x}{x_0}}\right)`$ (30)
$`{\displaystyle \frac{3\left(x_1s\right)^2}{x_1\left(x_1x_2\right)\left(x_1x_3\right)}}\mathrm{ln}\left({\displaystyle \frac{xx_1}{x_0x_1}}\right){\displaystyle \frac{3\left(x_2s\right)^2}{x_2\left(x_2x_1\right)\left(x_2x_3\right)}}\mathrm{ln}\left({\displaystyle \frac{xx_2}{x_0x_2}}\right)`$
$`{\displaystyle \frac{3\left(x_3s\right)^2}{x_3\left(x_3x_2\right)\left(x_3x_1\right)}}\mathrm{ln}\left({\displaystyle \frac{xx_3}{x_0x_3}}\right)],`$
where $`sa/M`$, $`x(r/M)^{1/2}`$, $`x_0(r_{ms}/M)^{1/2}`$, $`x_1`$, $`x_2`$, and $`x_3`$ are the three roots of $`x^33x+2s=0`$

$`x_1`$ $`=`$ $`2\mathrm{cos}\left({\displaystyle \frac{1}{3}}\mathrm{arccos}s{\displaystyle \frac{\pi }{3}}\right),`$ (31a)
$`x_2`$ $`=`$ $`2\mathrm{cos}\left({\displaystyle \frac{1}{3}}\mathrm{arccos}s+{\displaystyle \frac{\pi }{3}}\right),`$ (31b)
$`x_3`$ $`=`$ $`2\mathrm{cos}\left({\displaystyle \frac{1}{3}}\mathrm{arccos}s\right).`$ (31c)
It can be checked that $`f(r=r_{ms})=0`$, and $`f(rr_{ms})3M/2r^2`$.


format: text

text:
# 1 Introduction

## 1 Introduction

The experimental studies of nucleus-nucleus (A+A) collisions at the SPS and RHIC provide a rich information on hadron observables (multiplicities and momentum spectra). An extensive experimental program is motivated by a possibility to create a new state of matter – the quark gluon plasma (QGP) – in a laboratory. Any theoretical description of the QGP requires a macroscopic approach: statistical mechanics to describe the static properties of the matter and hydrodynamics to describe the system evolution. The equilibrium hadron gas (HG) model describes remarkably well the hadron multiplicities measured in Pb+Pb and Au+Au collisions at top SPS ($`\sqrt{s}=17`$ GeV per nucleon pair)  and RHIC ($`\sqrt{s}=130`$ GeV)  energies, where the creation of QGP is expected.

The HG model assumes the following formula for the hadron thermal multiplicities in the grand canonical ensemble (g.c.e.):

$$N_j=\frac{d_jV}{2\pi ^2}_0^{\mathrm{}}p^2𝑑p\left[\mathrm{exp}\left(\frac{\sqrt{p^2+m_j^2}\mu _j}{T}\right)\pm 1\right]^1,$$
(1)
where $`V`$ and $`T`$ are the system volume and temperature, respectively, $`m_j`$ and $`d_j`$ denote particle masses and degeneracy factors. The particle chemical potential $`\mu _j`$ in Eq.(1) is defined as $`\mu _j=b_j\mu _B+s_j\mu _S+c_j\mu _C`$ , where $`b_j,s_j,c_j`$ denote the baryonic number strangeness and charm of particle $`j`$. The baryonic chemical potential $`\mu _B`$ regulates the baryonic density of the HG system whereas strange $`\mu _S`$ and charm $`\mu _C`$ chemical potentials should be found from the requirement of zero value for the total strangeness and charm in the system (in our consideration we neglect small effects of a non-zero electrical chemical potential). The total multiplicities $`N_j^{tot}`$ in the HG model include the resonance decay contributions:

$$N_j^{tot}=N_j+\underset{R}{}Br(Rj)N_R,$$
(2)
where $`Br(Rj)`$ are the corresponding decay branching ratios. The hadron yield ratios $`N_j^{tot}/N_i^{tot}`$ in the g.c.e. are the functions of $`T`$ and $`\mu _B`$ and are independent of the volume parameter $`V`$.

The temperature parameters extracted from fitting the multiplicity data are approximately the same for the SPS and RHIC energies: $`T=170\pm 10`$ MeV . This value is close to an estimate of the temperature $`T_c`$ for the QGP–HG transition obtained in Lattice QCD simulations at zero baryonic density (see e.g. ). One may therefore argue that QGP formed at the early stage of A+A reaction hadronizes into a locally equilibrated HG, and the chemical composition of this hadron gas is weakly affected by the subsequent hadron rescatterings.

## 2  Statistical coalescence model of $`J/\psi `$ production

The production of charmonium states $`J/\psi `$ and $`\psi ^{}`$ have been measured in A+A collisions at CERN SPS over the last 15 years by the NA38 and NA50 Collaborations. These studies were mainly motivated by the theoretical suggestion  to use $`J/\psi `$ meson as a probe for deconfinement in A+A collisions.

Recently the thermal model  and the statistical coalescence model (SCM)  (see also ) for the charmonium production in A+A collisions were formulated. The total $`J/\psi `$ multiplicity, $`N_{J/\psi }^{tot}`$, in the thermal model  is given by Eq. (2), where $`N_{J/\psi }`$, $`N_\psi ^{}`$, $`N_{\chi _1}`$, $`N_{\chi _2}`$ are calculated according to Eq.(1) and $`Br(\psi ^{})0.54`$, $`Br(\chi _1)0.27`$, $`Br(\chi _2)0.14`$ are the decay branching ratios of the excited charmonium states into $`J/\psi `$. The thermal model  predicts that at high collision energies the $`J/\psi `$ to $`\pi `$ ratio is independent of $`\sqrt{s}`$ and the number of nucleon participants $`N_p`$. This is because both $`J/\psi `$ and $`\pi `$ multiplicities are proportional to the system volume and the hadronization temperature is expected to be approximately constant at high collision energies.

The SCM  assumes that charmonium states are formed at the hadronization stage. This is similar to the thermal model . However, in the SCM the charmonium states are produced via a coalescence of $`c`$ and $`\overline{c}`$ (anti)quarks created by the hard parton collisions at the early stage of A+A reaction. The number of created $`c\overline{c}`$ pairs, $`N_{c\overline{c}}`$, in hard parton collisions, differs in general from the result expected in the equilibrium HG. One needs then a charm enhancement factor $`\gamma _c`$ 111 This is formally analogous to the introduction of the strangeness suppression factor $`\gamma _s<1`$  in the HG model. to adjust the thermal HG results to the required average number, $`N_{c\overline{c}}`$, of $`c\overline{c}`$ pairs. The open charm hadron yield is enhanced by a factor $`\gamma _c`$ and charmonium yield by a factor $`\gamma _c^2`$ in comparison with the equilibrium HG predictions. This leads to a difference between the thermal model and SCM predictions for charmonia multiplicities.

The canonical ensemble (c.e.) formulation of the SCM is :

$$N_{c\overline{c}}=\frac{1}{2}\gamma _cN_O\frac{I_1(\gamma _cN_O)}{I_0(\gamma _cN_O)}+\gamma _c^2N_H,$$
(3)
where $`N_H`$ and $`N_O`$ are the HG multiplicities of all particles with hidden and open charm, respectively. The canonical suppression factor $`I_1/I_0`$ in Eq.(3) is due to the exact charm conservation (see e.g., , $`I_0`$ and $`I_1`$ are the modified Bessel functions).

If $`N_{c\overline{c}}`$ is known, Eq.(3) can be used to find the charm enhancement factor $`\gamma _c`$ and calculate then the $`J/\psi `$ multiplicity:

$$J/\psi =\gamma _c^2N_{J/\psi }^{tot},$$
(4)
where $`N_{J/\psi }^{tot}`$ is the total (thermal plus excited charmonium decays) HG $`J/\psi `$ multiplicity.

Note that the second term in the right-hand side of Eq.(3) gives only a tiny correction to the first term, i.e. most of the created $`c\overline{c}`$ pairs are transformed into the open charm hadrons. Eqs.(3,4) lead then to:

$$J/\psi N_{c\overline{c}}\frac{N_{J/\psi }^{tot}}{(N_O/2)^2},N_{c\overline{c}}<<1;$$
(5)
$$J/\psi N_{c\overline{c}}^2\frac{N_{J/\psi }^{tot}}{(N_O/2)^2},N_{c\overline{c}}>>1.$$
(6)
Eq.(3) assumes an exact conservation of $`N_cN_{\overline{c}}0`$ and ‘statistical fluctuations’ of $`N_{c\overline{c}}(N_c+N_{\overline{c}})/2`$ numbers. A more accurate treatment with ‘dynamical’ (Poisson-like) distribution of $`N_{c\overline{c}}`$ leads to :

$$J/\psi N_{c\overline{c}}\left(1+N_{c\overline{c}}\right)\frac{N_{J/\psi }^{tot}}{(N_O/2)^2}.$$
(7)
Eq.(7) coincides with Eqs.(3,4) at both $`N_{c\overline{c}}<<1`$ and $`N_{c\overline{c}}>>1`$ limits (see Eqs.(5,6)). This is because at $`N_{c\overline{c}}<<1`$ the probabilities to create zero, $`P(0)1N_{c\overline{c}}`$, and one, $`P(1)N_{c\overline{c}}`$, $`c\overline{c}`$ pairs are only important. On the other hand, at $`N_{c\overline{c}}>>1`$ the statistical and dynamical fluctuations of $`N_{c\overline{c}}`$ both obey the Poisson law distribution.

## 3  $`J/\psi `$ suppression in Pb+Pb collisions at the SPS

The centrality dependence of $`N_{c\overline{c}}`$ in A+B nucleus-nucleus collisions can be calculated in Glauber’s approach, $`N_{c\overline{c}}(b)=\sigma _{c\overline{c}}^{NN}T_{AB}(b)`$, where $`b`$ is the impact parameter, $`T_{AB}(b)`$ is the nuclear overlap function and $`\sigma _{c\overline{c}}^{NN}`$ is the $`c\overline{c}`$ production cross section for nucleon-nucleon collisions. As discussed in Ref., the deconfined medium can substantially modify charm production at the SPS, i.e., $`\sigma _{c\overline{c}}^{NN}`$ in A+B collisions can be different from the corresponding cross section measured in a nucleon-nucleon collision experiment. The present analysis treats $`\sigma _{c\overline{c}}^{NN}`$ at the SPS energy as a free parameter. Its value is fixed by fitting the NA50 data.

In the NA50 experiments  the Drell-Yan muon pair multiplicity (either measured or calculated from the minimum bias data) is used as a reference for the $`J/\psi `$ “suppression pattern”. The number of Drell-Yan pairs is also proportional to the number of primary nucleon-nucleon collisions: $`DY^{}(b)=\sigma _{DY^{}}^{NN}T_{AB}(b)`$, where $`\sigma _{DY^{}}^{NN}`$ is the nucleon-nucleon production cross section of $`\mu ^+\mu ^{}`$ Drell-Yan pairs. The prime means that the pairs should satisfy the kinematical conditions of the NA50 spectrometer. As the Drell-Yan cross section is isospin dependent, an average value is used: $`\sigma _{DY^{}}^{NN}=\sigma _{DY^{}}^{AB}/(AB)`$. For the case of Pb+Pb collisions, $`A=B=208`$ and $`\sigma _{DY^{}}^{PbPb}=1.49\pm 0.13`$ $`\mu `$b .

The quantity to be studied is the ratio

$`R_{DY}(b){\displaystyle \frac{\eta B_{\mu \mu }^{J/\psi }J/\psi (b)}{DY^{}(b)}}`$
$`=`$ $`\eta B_{\mu \mu }^{J/\psi }{\displaystyle \frac{\sigma _{c\overline{c}}^{NN}}{\sigma _{DY^{}}^{NN}}}{\displaystyle \frac{\left(1+\sigma _{c\overline{c}}^{NN}T_{AB}(b)\right)}{N_p(b)}}{\displaystyle \frac{n_{J/\psi }^{tot}(T,\mu _B)n_B(T,\mu _B)}{(n_O(T,\mu _B)/2)^2}}.`$
$`B_{\mu \mu }^{J/\psi }0.0588`$ is the decay probability of $`J/\psi `$ into $`\mu ^+\mu ^{}`$. Only the fraction $`\eta `$ of $`\mu ^+\mu ^{}`$ pairs satisfying the kinematical conditions of the NA50 spectrometer can be registered. We treat $`\eta `$ in Eq.(3) as one more model free parameter.

In the NA50 experiment, the neutral transversal energy $`E_T`$ of produced particles was used to measure centrality of the collisions. This variable, however, provides a reliable measure of the centrality only if it does not exceed a certain maximum value: $`E_T<100`$ GeV (see Ref.). The average number of participants is a linear function of the transversal energy, $`\overline{N}_p=E_T/q`$, in the domain $`E_T<100`$ GeV ($`q0.274`$ GeV). At $`E_T>100`$ GeV the number of $`\overline{N}_p`$ does not change essentially as $`E_T`$ grows. Therefore, the data at $`E_T>100`$ GeV do not represent a centrality dependence of the $`J/\psi `$ suppression pattern but rather its dependence on fluctuations of the stopping energy at approximately fixed number of participants. The influence of such fluctuations on $`J/\psi `$ multiplicity can be studied in the framework of the SCM model .

Fig. 1: The dependence of the $`J/\psi `$ to Drell-Yan ratio on the transversal energy. The normal nuclear suppression curve is obtained at $`\sigma _{abs}=6.4`$ mb. The SCM lines are calculated using Eq.(3). The vertical line shows the applicability domain of the SCM, $`N_p>100`$.

The SCM fit of the NA50 data with Eq.(3) is shown in Fig. 1 (see Refs.  for details). The model parameters are:

$$\sigma _{c\overline{c}}^{NN}35.7\mu \text{b},\eta 0.13,$$
(9)
with $`T`$ and $`\mu _B`$ fixed by hadron yield data. The model requires rather large enhancement for the open charm production (up to a factor of about $`3.5`$ within the rapidity window of the NA50 spectrometer). A direct measurement of the open charm in Pb+Pb at CERN SPS by NA60 would allow to test this prediction.

The SCM does not describe the NA50 data at small values of $`E_T`$ (the peripheral collisions with $`E<27`$ GeV). This can be also seen from the $`\psi ^{}`$ data. For $`T170`$ MeV the value of the thermal ratio is $`\psi ^{}/J/\psi 0.04`$. This is in agreement with data in Pb+Pb collisions at SPS for $`N_p>100`$, but is in contradiction with data in p+p, p+A and very peripheral A+A collisions. This fact was first noticed in Ref. .

## 4 $`J/\psi `$ enhancement in A+A Collisions at the RHIC

The number of directly produced $`c\overline{c}`$ pairs at the RHIC energies can be estimated in the pQCD approach and used then as the input for the SCM. The pQCD calculations for $`c\overline{c}`$ production cross sections were first done in Ref.. For the cross section $`\sigma (ppc\overline{c})`$ of the charm production in p+p collisions we use the results presented in Ref.. This leads to the value of $`\sigma (ppc\overline{c})0.35`$ mb at $`\sqrt{s}=200`$ GeV and the $`\sqrt{s}`$-dependence of the cross section for $`\sqrt{s}=10÷200`$ GeV is parameterized as :

$$\sigma (ppc\overline{c})=\sigma _0\left(1\frac{M_0}{\sqrt{s}}\right)^\alpha \left(\frac{\sqrt{s}}{M_0}\right)^\beta ,$$
(10)
with $`\sigma _03.392`$ $`\mu `$b, $`M_02.984`$ GeV, $`\alpha 8.185`$ and $`\beta 1.132`$.

The number of produced $`c\overline{c}`$ pairs in A+A collisions is proportional to the number of primary N+N collisions, $`N_{coll}^{AA}`$, which in turn is proportional to $`N_p^{4/3}`$ :

$$N_{c\overline{c}}=N_{coll}^{AA}(N_p)\frac{\sigma (ppc\overline{c})}{\sigma _{NN}^{inel}}C\sigma (ppc\overline{c})N_p^{4/3},$$
(11)
where $`\sigma _{NN}^{inel}30`$ mb is the inelastic N+N cross sections, $`C11`$ barn-1.

The results of the SCM can be studied analytically according to Eqs.(5,6) in the limiting cases of small and large numbers of $`N_{c\overline{c}}`$. For $`N_{c\overline{c}}<<1`$ one finds:

$$R\frac{J/\psi }{N_{c\overline{c}}}\frac{4N_{J/\psi }^{tot}}{N_O^2}\frac{1}{V}\frac{1}{\pi }\left(\sqrt{s}\right)^{1/2}N_p^1,$$
(12)
where we use the energy dependence of the pion multiplicity per nucleon participant $`\pi /N_p(\sqrt{s})^{1/2}`$  which approximately works in the SPS–RHIC energy region. The behavior (12) corresponds to the $`J/\psi `$ suppression: the ratio $`R`$ decreases with increasing of both $`\sqrt{s}`$ and $`N_p`$. This takes place at the SPS: this energy is still too “low” as $`N_{c\overline{c}}<1`$ even in the most central Pb+Pb collisions. However, the behavior of the $`J/\psi `$ to $`N_{c\overline{c}}`$ ratio is changed dramatically at the RHIC energies  (see Fig. 2 and Ref.  for details).

Fig. 2. Left: The energy dependence of the $`J/\psi `$ to $`N_{c\overline{c}}`$ ratio in central Au+Au collisions. Points are the predictions of the SCM for the RHIC energies: $`\sqrt{s}=56,130,200`$ GeV. Right: The $`N_p`$-dependence of the $`J/\psi `$ to $`N_{c\overline{c}}`$ ratio. The lines are the predictions of the SCM. The dashed line corresponds to $`\sqrt{s}=56`$ GeV (the $`J/\psi `$ suppression), the solid line corresponds to $`\sqrt{s}=200`$ GeV (the $`J/\psi `$ enhancement).

In central Au+Au collisions at $`\sqrt{s}=200`$ GeV the expected value of $`N_{c\overline{c}}`$ is essentially larger than unit. For $`N_{c\overline{c}}>>1`$ one finds ($`\beta 1.1`$):

$$R\frac{J/\psi }{N_{c\overline{c}}}\frac{2N_{c\overline{c}}}{N_O}\frac{N_{c\overline{c}}}{V}\frac{N_{c\overline{c}}}{\pi }(\sqrt{s})^{\beta 1/2}N_p^{1/3}.$$
(13)
## 5 The $`m_T`$-spectra of $`J/\psi `$ and $`\psi ^{}`$ and QGP hadronization

In Ref. we formulated the hypothesis that the kinetic freeze-out of $`J/\psi `$ and $`\psi ^{}`$ mesons takes place directly at hadronization. The effect of rescattering in the hadronic phase was recently studied within a “hydro + cascade” approach . A+A collisions are considered there to proceed in three stages: hydrodynamic QGP expansion (“hydro”), transition from QGP to HG (“switching”) and the stage of hadronic rescatterings and resonance decays (“cascade”). The switching from hydro to cascade takes place at $`T=T_c`$, where the spectrum of hadrons leaving the surface of the QGP–HG transition is taken as an input for the subsequent cascade calculations. The results of Refs.  suggest that the transverse momentum spectra of $`\mathrm{\Omega }`$ is only weakly affected during the cascade stage 222The corresponding calculation for charmonia are not yet performed within this model. and give therefore a straightforward measure of the collective hydro velocity at the ‘switching’ surface $`T=T_c`$. In Ref.  we demonstrated that in Pb+Pb collisions at 158 A$``$GeV the $`m_T`$-spectra of $`\mathrm{\Omega }^\pm `$  can be explained simultaneously with the $`m_T`$-spectra of $`J/\psi `$ and $`\psi ^{}`$ mesons  using the same set of the hadronization parameters.

Assuming the fluid freeze-out at constant temperature $`T`$, the transverse mass spectrum of $`i`$-th hadron species in cylindrically symmetric and longitudinally boost invariant fluid expansion equals approximately to (see Ref.  for further references and details):

$$\frac{dN_i}{m_Tdm_T}\sqrt{m_T}\mathrm{exp}\left(\frac{m_T(1+\frac{1}{2}\overline{v}_T^2)}{T}\right)I_0\left(\frac{p_T\overline{v}_T}{T}\right),$$
(14)
where $`m_T=(m_i^2+p_T^2)^{1/2}`$ and $`\overline{v}_T`$ is the average transverse flow velocity.

Fig. 3: The $`m_T`$-spectra measured at midrapidity in Pb+Pb at 158 A$``$GeV by WA97  for $`\mathrm{\Omega }`$ and by NA50  for $`J/\psi `$ and $`\psi ^{}`$ are presented (in arbitrary units) versus $`m_Tm`$. The solid lines correspond to Eq.(14) with $`T=170`$ MeV and $`\overline{v}_T=0.19`$.

In Fig. 3 we present the fit with Eq.(14) of the measured $`m_T`$-spectra. The temperature is fixed as $`T=170`$ MeV and $`\overline{v}_T`$ is considered as a free parameter. The shapes of all spectra are simultaneously reproduced at $`\overline{v}_T=0.19\pm 0.02`$.

## 6  Conclusions

Statistical hadronization of the QGP is probably an important source of $`J/\psi `$ production . Within the SCM  the NA50 data on the $`J/\psi `$ production in Pb+Pb at 158 A$``$GeV can be fitted  (see Fig. 1) for central collisions $`N_p>100`$. A large enhancement of the open charm over an extrapolation from the p+p data is however required . A direct measurement of the open charm in Pb+Pb at CERN SPS by NA60 would allow to test the SCM selfconsistency.

The SCM predicts that the $`J/\psi `$ suppression at the SPS should be changed into the $`J/\psi `$ enhancement in central Au+Au collisions at the RHIC energies  (see Fig. 2).

The shapes of $`m_T`$-spectra for $`\mathrm{\Omega }`$, $`J/\psi `$ and $`\psi ^{}`$ are simultaneously reproduced  (see Fig. 3) in the hydrodynamical picture of the QGP hadronization with $`T=170`$ MeV, $`\overline{v}_T=0.19\pm 0.02`$. This supports the hypothesis that formation and the kinetic freeze–out of charmonia occurs at the hadronization .

Acknowledgments. I am thankful to K.A. Bugaev, M. Gaździcki, W. Greiner, A.P. Kostyuk, L. McLerran and H. Stöcker for fruitful collaboration. I am also thankful to F. Becattini, L. Bravina, P. Braun-Munzinger, J. Cleymans, A. Dumitru, L. Gerland, D. Kharzeev, I.N. Mishustin, G.C. Nayak, K. Redlich, Yu.M. Sinyukov, J. Stachel, D. Teaney and Nu Xu for comments and discussions. The financial support from the Humboldt Foundation is acknowledged. The research described in this publication was made possible in part by Award # UP1-2119 of the U.S. Civilian Research and Development Foundation for the Independent States of the Former Soviet Union (CRDF) and INTAS grant 00-00366.


text:
# Strong Coupling Limit of the Kardar-Parisi-Zhang Equation in 2+1 Dimensions

## I Introduction

Because of technical importance and fundamental interest, a great deal of efforts have been devoted to the understanding of the mechanism of thin-film growth and the kinetic roughening of growing surfaces in various growth techniques. Analytical and numerical treatments of simple growth models suggest that, quite generally, the height fluctuations have a self-similar character and their average correlations exhibit a dynamical scaling form. Numerous theoretical models have been proposed, of which the simplest nontrivial example is the Kardar-Parisi-Zhang equation

$$\frac{h}{t}\frac{\alpha }{2}(h)^2=\nu ^2h+f(x,y,t),$$
(1)
where $`h(𝐱,𝐭)`$ specifies the height of the surface at point $`𝐱`$ and $`f`$ is a zero-mean, statistically homogeneous, white in time Gaussian process with covariance $`f(𝐱,t)f(𝐱^{},t^{})=2D_0\delta (tt^{})D(𝐱𝐱^{})`$. Typically the spatial correlation of forcing is considered to be as a delta function, mimicking the short range correlation. Here we consider the spatial correlation as $`D(𝐱𝐱^{})=\frac{\mathrm{𝟏}}{\pi \sigma ^\mathrm{𝟐}}\mathrm{exp}(\frac{(𝐱𝐱^{})^\mathrm{𝟐}}{\sigma ^\mathrm{𝟐}}),`$ where $`\sigma `$ is the variance of $`D(𝐱𝐱^{})`$. When the variance $`\sigma `$ is much less than the system size $`L`$ , i.e. $`\sigma L`$, the model represents a short range character for the forcing. This same equation is believed to describe the statistics of directed polymers in a random medium , where $`h`$ has the meaning of the free energy. The average force on the interface is not essential and can be removed by a simple shift in $`h`$. Every term in the eq.(1) involves a specific physical phenomenon contributing to the surface evolution. The parameters $`\nu `$, $`\alpha `$ and $`D_0`$ ( and $`\sigma `$) are related to surface relaxation, lateral growth and the noise strength, respectively. The effective coupling constant for the KPZ equation is given by $`g=\frac{2\alpha ^2D_0}{\nu ^3}`$. The limit $`g\mathrm{}`$ is known as the strong coupling limit. Despite the intense effort in recent years, the properties of the strong coupling phase are rather poorly understood. Only the critical exponents of the strong-coupling regime ($`g\mathrm{}`$ or $`\nu 0`$) are known in 1+1 dimensions and their values in higher dimensions as well as properties of the roughening transition have been known only numerically and by the various approximative schemes.
Theoretical richness of the KPZ model is partly due to its close relationships with other areas of statistical physics. It is shown that there is a mapping between the equilibrium statistical mechanics of a two dimensional smectic-A liquid crystal onto the non-equilibrium dynamics of the (1+1)- dimensional stochastic KPZ equation . It has been shown  that, one can map the kinetics of the annihilation process $`A+B0`$ with driven diffusion onto the (1+1)-dimensional KPZ equation. Also the KPZ equation is closely related to the dynamics of a Sine-Gordon chain , the driven-diffusion equations , and directed paths in the random media  and charge density waves , dislocations in disordered solids , formation of large-scale structure in the universe , Burgers turbulence  and etc.

In this paper we attempt to clarify the short time physics of KPZ equation in the strong coupling limit by presenting an exact solution of the model in 2+1 dimensions. It turns out that eq.(1) exhibits a finite time singularity at time scale $`t_c`$. Using the master equation method it is shown that $`t_c`$ can be expressed in terms of the characteristics of forcing as, $`t_c=\frac{1}{4}(\frac{\pi }{8\alpha ^2D_0})^{1/3}\sigma ^2`$. Also we derive the time dependence of the moments $`<(h\overline{h})^n>`$ for time scales before creation of singularities.

The main properties of the KPZ equation in the strong-coupling limit are as follow: (i) In the limit of $`\nu 0`$ the unforced KPZ (Burgers) equation develops singularities for the given dimension. In one spatial dimension it develops sharp valleys (point shocks). In the singular points the height gradients are not continuous . In two spatial dimensions KPZ equation develops three types of singularities, the first singularities are finite sharp valley lines (shock lines) across which the height gradients are discontinuous. The second type is the end point of the sharp valley lines which separates the regular points and singular region and is called a kurtoparabolic point. As time increases these sharp valley lines hit each other and crossing point of two valley lines (shock lines) produces a valley node (shock node). Generically kurtoparabolic points disappear at large times and only a network of sharp valley lines survive . A complete classification of the singularities of KPZ (Burgers) equation in two and three dimensions, by considering the metamorphoses of singularities as time elapses, has been done in . (ii) For white in time and smooth in space forcing, the sharp valley lines (cusp lines ) are smooth curves and the singularity structures of the forced case is similar to unforced problem. It is shown in  that in the stationary state the sharp valley lines produce a curvilinear hexagon lattice in two dimensions and, therefore, one finds a hexagonal tiling of singularity lines (valley lines).

Here we shall adapt the master equation approach which enables us to investigate the KPZ equation in the strong coupling limit . It is one of the advantages of this method that allows us to write all the nonlinearities due to the nonlinear term $`\frac{\alpha }{2}(h)^2`$ in a closed form. When $`\sigma `$ is finite, the nonlinear term develops finite time singularities for the height derivatives in the limit ($`\nu 0`$). Therefore, one would distinguish between different time regimes before and after the sharp valley formation. Starting from a flat initial condition, i.e. $`h(x,0)=0,h_x=u(x,0)=0`$ and $`h_y=v(x,0)=0`$, we aim at tracing the time evolution given by the KPZ equation with $`\nu 0`$. Only for time scales before the sharp valley formation ($`t<t_c`$) the limit $`\nu 0`$ and $`\nu =0`$ are equivalent. For $`t>t_c`$, the anomalous contributions of the diffusion terms to the equation of probability distribution function will break this equivalence . Therefore, for $`t<t_c`$, one can impose $`\nu =0`$ in the dynamics and try to find realisable statistical solutions for the problem. We show that it is possible to find finite time realisable solutions for the probability density function of the height field at least in the early steps of the evolution of the surface. We attribute the finite life time of the realisable ensemble of statistical solutions to the fact that after a finite time the derivatives of the field $`h(x,t)`$ become singular and , therefore, the diffusion term contributions are no more negligible.

## II calculation of the height moments before the singularity formation

Let us consider the Kardar-Parisi-Zhang equation in 2+1 dimensions,

$$h_t(x,y,t)\frac{\alpha }{2}(h_x^2+h_y^2)=\nu (h_{xx}+h_{yy})+f(x,y,t).$$
(2)
Defining $`h_x(x,y,t)=u(x,y,t)`$ and $`h_y(x,y,t)=v(x,y,t)`$, differentiating the above equation with respect to $`x`$ and $`y`$ and neglecting the viscosity term in the limit of $`\nu 0`$ before the creation of singularities we get

$`u_t(x,y,t)`$ $`=`$ $`\alpha (u(x,y,t)u_x(x,y,t)+v(x,y,t)u_y(x,y,t))+f_x(x,y,t)`$ (3)
$`v_t(x,y,t)`$ $`=`$ $`\alpha (u(x,y,t)v_x(x,y,t)+v(x,y,t)v_y(x,y,t))+f_y(x,y,t).`$ (4)
Now introducing $`\mathrm{\Theta }`$ as

$`\mathrm{\Theta }=`$ (5)
$`\mathrm{exp}\left(i\lambda (h(x,y,t)\overline{h}(t))i\mu _1u(x,y,t)i\mu _2v(x,y,t)\right).`$ (6)
enables us to express the generating function as $`Z(\lambda ,\mu _1,\mu _2,x,y,t)=\mathrm{\Theta }`$. Using the KPZ equation and its differentiations with respect to $`x`$ and $`y`$ and neglecting the viscosity term, in the limit of $`\nu 0`$ before the creation of singularities, we get the following expression for the time evolution of $`Z(\lambda ,\mu _1,\mu _2,x,y,t)`$

$`Z_t`$ $`=`$ $`i\gamma (t)\lambda Zi\lambda {\displaystyle \frac{\alpha }{2}}u^2\mathrm{\Theta }i\lambda {\displaystyle \frac{\alpha }{2}}v^2\mathrm{\Theta }`$ (7)
$``$ $`i\alpha \mu _1uu_x\mathrm{\Theta }i\alpha \mu _1vv_x\mathrm{\Theta }i\alpha \mu _2uu_y\mathrm{\Theta }i\alpha \mu _2vv_y\mathrm{\Theta }`$ (8)
$``$ $`i\lambda f(x,y,t)\mathrm{\Theta }i\mu _1f_x(x,y,t)\mathrm{\Theta }i\mu _2f_y(x,y,t)\mathrm{\Theta }`$ (9)
where $`\gamma (t)=\overline{h}_t=\frac{\alpha }{2}<u^2+v^2>,k(𝐱𝐱^{})=2D_0D(𝐱𝐱^{}),k^{\prime \prime }(0,0)=k_{xx}(0,0)=k_{yy}(0,0)`$. To simplify the above expression, by taking the derivative of $`Z`$ with respect to $`x`$ and $`y`$, and using the homogeneity conditions, we obtain

$`Z_x`$ $`=`$ $`(i\lambda ui\mu _1u_xi\mu _2u_y)\mathrm{\Theta }=0`$ (10)
$`Z_y`$ $`=`$ $`(i\lambda vi\mu _1v_xi\mu _2v_y)\mathrm{\Theta }=0,`$ (11)
from which it can be easily derived

$$i\frac{}{\mu _1}(i\mu _1u_xi\mu _2u_y)\mathrm{\Theta }=u_x\mathrm{\Theta }i\mu _1uu_x\mathrm{\Theta }i\mu _2uu_y\mathrm{\Theta }.$$
(12)
From the eqs. (10) and (12) we have

$$u_x\mathrm{\Theta }i\mu _1uu_x\mathrm{\Theta }i\mu _2uu_y\mathrm{\Theta }=\lambda \frac{}{\mu _1}u\mathrm{\Theta }=i\lambda Z_{\mu _1\mu _1}$$
(13)
from which we find the following expression for $`i\alpha \mu _1uu_x\mathrm{\Theta }i\alpha \mu _2uu_y\mathrm{\Theta }`$

$$i\alpha \mu _1uu_x\mathrm{\Theta }i\alpha \mu _2uu_y\mathrm{\Theta }=i\alpha \lambda Z_{\mu _1\mu _1}\alpha u_x\mathrm{\Theta }.$$
(14)
Similarly for $`i\alpha \mu _1vv_x\mathrm{\Theta }i\alpha \mu _2vv_y\mathrm{\Theta }`$ we find

$$i\alpha \mu _1vv_x\mathrm{\Theta }i\alpha \mu _2vv_y\mathrm{\Theta }=i\alpha \lambda Z_{\mu _2\mu _2}\alpha v_y\mathrm{\Theta }$$
(15)
Now from eqs. (14) and (15) and by using the Novikov’s theorem the eq.(7) gives

$`Z_t`$ $`=`$ $`i\gamma (t)\lambda Z+i\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _1\mu _1}+i\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _2\mu _2}i\alpha \lambda Z_{\mu _1\mu _1}i\alpha \lambda Z_{\mu _2\mu _2}`$ (16)
$``$ $`\alpha u_x\mathrm{\Theta }\alpha v_y\mathrm{\Theta }\lambda ^2k(0,0)Z+\mu _1^2k^{\prime \prime }(0,0)Z+\mu _2k^{\prime \prime }(0,0)Z.`$ (17)
The term $`u_x\mathrm{\Theta }`$ can be evaluated as follows

$`u_x\mathrm{\Theta }`$ $`=`$ $`{\displaystyle \frac{i}{\mu _1}}\mathrm{\Theta }_x+{\displaystyle \frac{i}{\mu _1}}(i\lambda u+i\mu _2v_x)\mathrm{\Theta }`$ (18)
$`=`$ $`i{\displaystyle \frac{\lambda }{\mu _1}}Z_{\mu _1}{\displaystyle \frac{\mu _2}{\mu _1}}v_x\mathrm{\Theta },`$ (19)
also for $`v_y\mathrm{\Theta }`$ we have

$$v_y\mathrm{\Theta }=i\frac{\lambda }{\mu _2}Z_{\mu _2}\frac{\mu _1}{\mu _2}u_y\mathrm{\Theta }.$$
(20)
In the appendix it is proved that $`u_y\mathrm{\Theta }=v_x\mathrm{\Theta }=0`$, so the equation governing the time evolution of $`Z`$ will become

$`Z_t`$ $`=`$ $`i\gamma (t)\lambda Zi\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _1\mu _1}i\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _2\mu _2}`$ (21)
$`+`$ $`i\alpha {\displaystyle \frac{\lambda }{\mu _1}}Z_{\mu _1}i\alpha {\displaystyle \frac{\lambda }{\mu _2}}Z_{\mu _2}`$ (22)
$``$ $`\lambda ^2k(0,0)Z+\mu _1^2k^{\prime \prime }(0,0)Z+\mu _2k^{\prime \prime }(0,0)Z.`$ (23)
In what follows we solve this partial differential equation by using flat initial condition, $`h(x,y,0)=u(x,y,0)=v(x,y,0)=0`$,which means

$$P(\stackrel{~}{h},u,v,0)=\delta (\stackrel{~}{h})\delta (u)\delta (v),$$
(24)
so that the initial condition for the generating function is $`Z(\lambda ,\mu _1,\mu _2,0)=1`$.
The solution of the eq.(21) can be factorized as

$$Z(\lambda ,\mu _1,\mu _2,t)=F_1(\lambda ,\mu _1,t)F_2(\lambda ,\mu _2,t)\mathrm{exp}(\lambda ^2k(0,0)t),$$
(25)
which by inserting eq.(25) in eq.(21) the following equation is obtained

$`F_{1}^{}{}_{t}{}^{}F_2+F_1F_{2}^{}{}_{t}{}^{}`$ $`=`$ $`i\gamma (t)\lambda F_1F_2i\lambda {\displaystyle \frac{\alpha }{2}}F_2F_{1}^{}{}_{\mu _1\mu _1}{}^{}i\lambda {\displaystyle \frac{\alpha }{2}}F_1F_{2}^{}{}_{\mu _2\mu _2}{}^{}`$ (26)
$`+`$ $`i\alpha {\displaystyle \frac{\lambda }{\mu _1}}F_2F_{1}^{}{}_{\mu _1}{}^{}i\alpha {\displaystyle \frac{\lambda }{\mu _2}}F_1F_{2}^{}{}_{\mu _2}{}^{}`$ (27)
$``$ $`\lambda ^2k(0,0)F_1F_2+\mu _1^2k^{\prime \prime }(0,0)F_1F_2+\mu _2k^{\prime \prime }(0,0)F_1F_2.`$ (28)
On the other hand, $`\gamma (t)=\overline{h}_t=\frac{\alpha }{2}u^2+v^2`$ and $`\frac{\alpha }{2}u^2=\frac{\alpha }{2}v^2=\alpha k^{\prime \prime }(0,0)t`$, so eq.(26) can be separated in terms of $`\mu _1`$ and $`\mu _2`$ and it can be seen easily that $`F_1(\lambda ,\mu _1,t)`$ and $`F_2(\lambda ,\mu _2,t)`$ are satisfied by the following 1+1-dimensional equation

$$F_t(\lambda ,\mu ,t)=i\lambda \frac{\alpha }{2}F_{\mu \mu }+i\alpha \frac{\lambda }{\mu }F_\mu +[\mu ^2k^{\prime \prime }(0,0)i\alpha \lambda k^{\prime \prime }(0,0)t]F,$$
(29)
which can be solved exactly with the initial condition $`F(\lambda ,\mu ,0)=1`$. The solution of eq.(29) is

$`F(\mu ,\lambda ,t)=(1\mathrm{tanh}^2(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t))`$ (32)
$`\mathrm{exp}[{\displaystyle \frac{5}{8}}\mathrm{ln}(1\mathrm{tanh}^4(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t))`$
$`+`$ $`{\displaystyle \frac{5}{4}}\mathrm{tanh}^1(\mathrm{tanh}^2(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t))\lambda ^2k(0,0)t`$ (34)
$``$ $`{\displaystyle \frac{1}{16}}\mathrm{ln}^2({\displaystyle \frac{1\mathrm{tanh}(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t)}{1+\mathrm{tanh}(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t)}})`$ (36)
$``$ $`{\displaystyle \frac{1}{2}}i\mu ^2\sqrt{{\displaystyle \frac{2ik_{xx}(0,0)}{\alpha \lambda }}}\mathrm{tanh}(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t)].`$ (38)
So, the solution of eq.(26) is

$`Z(\lambda ,\mu _1,\mu _2,t)=F(\lambda ,\mu _1,t)F(\lambda ,\mu _2,t)\mathrm{exp}(\lambda ^2k(0,0)t)`$ (39)
$`=(1\mathrm{tanh}^2(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t))^2\mathrm{exp}[{\displaystyle \frac{5}{4}}\mathrm{ln}(1\mathrm{tanh}^4(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t))`$ (40)
(41)
$`+{\displaystyle \frac{5}{2}}\mathrm{tanh}^1(\mathrm{tanh}^2(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t))\lambda ^2k(0,0)t`$ (42)
(43)
$`{\displaystyle \frac{1}{8}}\mathrm{ln}^2({\displaystyle \frac{1\mathrm{tanh}(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t)}{1+\mathrm{tanh}(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t)}})`$ (44)
(45)
$`{\displaystyle \frac{1}{2}}i(\mu _1^2+\mu _2^2)\sqrt{{\displaystyle \frac{2ik_{xx}(0,0)}{\alpha \lambda }}}\mathrm{tanh}(\sqrt{2ik_{xx}(0,0)\alpha \lambda }t)].`$ (46)
Using the generating function, the probability distribution function (PDF) of height fluctuation can be derived by inverse Fourier transformation as

$`P(\stackrel{~}{h},u,v,t)=`$ (47)
$`{\displaystyle \frac{d\lambda }{2\pi }\frac{d\mu _1}{2\pi }\frac{d\mu _2}{2\pi }\mathrm{exp}\left(i\lambda \stackrel{~}{h}+i\mu _1u+i\mu _2v\right)Z(\lambda ,\mu _1,\mu _2,t)}.`$ (48)
By expanding the solution of the generating function in powers of $`\lambda `$, all the $`(h\overline{h})^n`$ moments can be derived. For instance, the first five moments before the sharp valley formation are

$`\stackrel{~}{h}^2`$ $`=`$ $`\left({\displaystyle \frac{k^2(0,0)}{\alpha k^{\prime \prime }(0,0)}}\right)^{2/3}\left[{\displaystyle \frac{2}{3}}({\displaystyle \frac{t}{t^{}}})^4+2{\displaystyle \frac{t}{t^{}}}\right]`$ (49)
$`\stackrel{~}{h}^3`$ $`=`$ $`{\displaystyle \frac{48}{45}}\left({\displaystyle \frac{k^2(0,0)}{\alpha k^{\prime \prime }(0,0)}}\right)({\displaystyle \frac{t}{t^{}}})^6`$ (51)
$`\stackrel{~}{h}^4`$ $`=`$ $`\left({\displaystyle \frac{k^2(0,0)}{\alpha k^{\prime \prime }(0,0)}}\right)^{4/3}\left[{\displaystyle \frac{44}{35}}({\displaystyle \frac{t}{t^{}}})^88({\displaystyle \frac{t}{t^{}}})^5+12({\displaystyle \frac{t}{t^{}}})^2\right]`$ (53)
$`\stackrel{~}{h}^5`$ $`=`$ $`\left({\displaystyle \frac{k^2(0,0)}{\alpha k^{\prime \prime }(0,0)}}\right)^{5/3}\left[{\displaystyle \frac{1216}{945}}({\displaystyle \frac{t}{t^{}}})^{10}+{\displaystyle \frac{64}{3}}({\displaystyle \frac{t}{t^{}}})^7\right]`$ (55)
where $`t_{}=(\frac{k(0,0)}{\alpha ^2k^{\prime \prime 2}(0,0)})^{1/3}`$. An important information content of the derived exact form is that they determine the time scale of the singularity formations. One should first check the realisability condition, i.e. $`P(h\overline{h},t)0`$. In fact the moment relations above indicate that different even order moments will get negative in some distinct characteristic time scales. A closer inspection into the even moment relations reveals that the higher the moments are, the smaller their characteristic time scales become. So the rate of decreasing tends to $`t_c=\frac{1}{4}t_{}`$ for very large even moments asymptotically, where $`t_{}=(\frac{k(0,0)}{\alpha ^2k^{\prime \prime 2}(0,0)})^{1/3}=(\frac{\pi }{8\alpha ^2D_0})^{1/3}\sigma ^2`$. Therefore, it is concluded that at this time the right tail of the probability distribution function (PDF) starts to become negative, which is the reminiscent of singularity creation. Physically, the connection between the negativity of the right tail and the singularity formation time scale is related to a multi-valued solution of KPZ (Burgers) equation for $`tt_c`$. It should be notified that eq.(16) preserves the normalisability, i.e. $`Z(0,0,0,t)=1`$, equivalent to $`_{\mathrm{}}^+\mathrm{}P(h\overline{h},u,v;t)d(h\overline{h})𝑑u𝑑v=1`$ for every time $`t`$. So the PDF of $`h\overline{h}`$ and its derivatives are always normalizable to unity. In the $`\nu =0`$ limit, for time scales greater than $`t_c`$ , the height field will become multi-valued in the valleys (the shock region), which is related to the $`P(h\overline{h})`$ left tail . The multiplicity of the height field will increase the probability measure in the PDF. Therefore, to compensate the exceeded measure related to the multi-valued solutions the right tail of the PDF tails should become negative. The singularities in the limit $`\nu 0`$ can be constructed from multi-valued solutions of the KPZ equation with $`\nu =0`$ by Maxwell cutting rule , which makes the discontinuity in the derivative of the height filed. When $`t>t_c`$, the contribution of the relaxation term in the limit of vanishing diffusion coefficient should also be considered in order to find a realisable probability density function of height field. In other words, the disregardance of the PDF equation’s diffusion term is only valid up to the time scales in which the singularities have not been developed yet.

Taking into account $`\alpha >0`$ and $`k^{\prime \prime }(0,0)<0`$, the odd order moments will be positive in time scales $`t<t_c`$ and the probability density function $`P(h\overline{h},t)`$ is positively skewed. Therefore, the probability distribution functions of height field will have a non zero skewness through its time evolution at least before the singularity formations time scale.The inverse Fourier transform of the generating function has been performed numerically to obtain the PDF $`P(h\overline{h})`$ form. To demonstrate the time scale of the singularity formation the PDF time evolution has been sketched numerically, Fig(1). As the system evolves in time, the formation of the first singularities leads to the right tail negativity in the PDF.

To summarise, we analyse the fully nonlinear KPZ equation in 2$`+`$1 dimensions forced with a Gaussian noise which is white in time and short range correlated in space. In the non-stationary regime when the sharp valley structures are not yet developed we find an exact form for the generating function of the joint fluctuations of height and height gradient. We determine the time scale of the sharp valley formation and the exact functional form of the time dependence in the height difference moments at any given order. In this paper, we have used the initial condition $`h(x,0)=0`$ and $`u(x,0)=0`$. This method enables us to determine the finite time scale of singularity formation in terms of statistical properties of the initial conditions. This calculation will be dealt with elsewhere. At this level we would remark on the stationary state of the KPZ equation in the strong coupling limit. Using the numerical results recently, it is shown that convergence to the statistical steady state is reached after a few turnover times . Therefore, in the stationary state the singularities will be fully developed and, therefore, at large $`t`$ or $`t\mathrm{}`$ we should take into account the relaxation contribution in the limit $`\nu 0`$. Contribution of the relaxation term in the rhs of the eq.(6) is,

$$\underset{\nu 0}{lim}\{i\nu \mu _1^2u\mathrm{\Theta }i\nu \mu _2^2v\mathrm{\Theta }\}.$$
(56)
In  we have calculated the finite contribution of this term in 1+1 dimensions and showed that it plays an important role in statistical theory of KPZ equation in the stationary state. It is also shown that all of the amplitudes of the moments $`<(h\overline{h})^n>`$ can be expressed in term of the relaxation term in the limit $`\nu 0`$. In 2+1 dimension calculation of finite contribution of the relaxation term is more complex ( because of complex structure of the singularities) and it is left for a future work. We believe that the analysis followed in that paper is also suitable for the zero temperature limit in the problem of directed polymer in the random potential with short range correlations . In the same direction the present method is applicable to the strong coupling regime of KPZ equation in higher dimensions ($`d>2`$) which is definitely an important step.
Acknowledgement
We thank J. Davoudi, F. Arash, M. Khorrami, A. Reza-Khani and A. Tabei for useful discussions. F.S and A.A.M. were partially supported by the Institute for Studies in Physics and Mathematics (IPM).

## III APPENDIX

In this appendix we present another way of deriving the height moments by introducing the second spatial derivatives of height field. Defining $`\mathrm{\Theta }(\lambda ,\mu _1,\mu _2,\eta _1,\eta _2,\eta _3,x,y,t)`$ as

$`\mathrm{\Theta }(\lambda ,\mu _1,\mu _2,\eta _1,\eta _2,\eta _3,x,y,t)=`$ $`\mathrm{exp}\{i\lambda (h(x,y,t)\overline{h}(t))i\mu _1u(x,y,t)i\mu _2v(x,y,t)`$ (58)
$`i\eta _1w(x,y,t)i\eta _2s(x,y,t)i\eta _3q(x,y,t)\},`$
the generating function is defined as

$$Z(\lambda ,\mu _1,\mu _2,\eta _1,\eta _2,\eta _3,x,y,t)=\mathrm{\Theta },$$
(59)
where $`u(x,y,t)=h_x(x,y,t)`$,$`v(x,y,t)=h_y(x,y,t)`$,$`w(x,y,t)=h_{xx}(x,y,t)`$, $`s(x,y,t)=h_{xy}(x,y,t)`$ and $`q(x,y,t)=h_{yy}(x,y,t)`$.

The evolution of $`h(x,y,t),u(x,y,t),v(x,y,t),w(x,y,t),s(x,y,t),q(x,y,t)`$ is given by the following equations

$`h_t`$ $`=`$ $`{\displaystyle \frac{\alpha }{2}}(u^2+v^2)+f(x,y,t)`$ (60)
$`u_t`$ $`=`$ $`\alpha (uw+vs)+f_x(x,y,t)`$ (61)
$`v_t`$ $`=`$ $`\alpha (us+vq)+f_y(x,y,t)`$ (62)
$`w_t`$ $`=`$ $`\alpha (w^2+s^2+uw_x+vs_x)+f_{xx}(x,y,t)`$ (63)
$`s_t`$ $`=`$ $`\alpha (ws+qs+vs_yuw_y)+f_{xy}(x,y,t)`$ (64)
$`q_t`$ $`=`$ $`\alpha (s^2+q^2+us_y+vq_y)+f_{yy}.`$ (65)
It follows from the above equations that the generating function $`Z`$ is the solution of the following equation

$`Z_t`$ $`=`$ $`i\gamma (t)\lambda Zi\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _1\mu _1}i\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _2\mu _2}+i\alpha Z_{\mu _1x}+i\alpha Z_{\mu _2y}`$ (66)
$``$ $`i\alpha Z_{\eta _1}i\alpha Z_{\eta _3}+i\alpha \eta _1Z_{\eta _1\eta _1}+i\alpha \eta _1Z_{\eta _2\eta _2}`$ (67)
$`+`$ $`i\alpha \eta _2Z_{\eta _1\eta _2}+i\alpha \eta _2Z_{\eta _2\eta _3}+i\alpha \eta _3Z_{\eta _2\eta _2}+i\alpha \eta _3Z_{\eta _3\eta _3}`$ (68)
$``$ $`i\lambda f(x,y,t)\mathrm{\Theta }i\mu _1f_x(x,y,t)\mathrm{\Theta }i\mu _2f_y(x,y,t)\mathrm{\Theta }`$ (69)
$``$ $`i\eta _1f_{xx}(x,y,t)\mathrm{\Theta }i\eta _2f_{xy}(x,y,t)\mathrm{\Theta }i\eta _3f_{yy}(x,y,t)\mathrm{\Theta },`$ (70)
in which $`\gamma (t)`$ is defined as $`\gamma (t)=\overline{h}_t`$ and the following identities have been used

$`\eta _1w_x\mathrm{\Theta }+\eta _2s_x\mathrm{\Theta }+\eta _3q_x\mathrm{\Theta }`$ $`=`$ $`iZ_xi\lambda Z_{\mu _1}i\mu _1Z_{\eta _1}i\mu _2Z_{\eta _2}`$ (71)
$`\eta _1w_y\mathrm{\Theta }+\eta _2s_y\mathrm{\Theta }+\eta _3q_y\mathrm{\Theta }`$ $`=`$ $`iZ_yi\lambda Z_{\mu _2}i\mu _1Z_{\eta _2}i\mu _2Z_{\eta _3}.`$ (72)
Now, using Novikov’s theorem we find

$`f(x,y,t)\mathrm{\Theta }`$ $`=`$ $`i\lambda k(0,0)Zi\eta _1k_{xx}(0,0)Zi\eta _3k_{xx}(0,0)Z`$ (73)
$`f_x(x,y,t)\mathrm{\Theta }`$ $`=`$ $`i\mu _1k_{xx}(0,0)Z`$ (74)
$`f_y(x,y,t)\mathrm{\Theta }`$ $`=`$ $`i\mu _2k_{xx}(0,0)Z`$ (75)
$`f_{xx}(x,y,t)\mathrm{\Theta }`$ $`=`$ $`i\lambda k_{xx}(0,0)Zi\eta _1k_{xxxx}(0,0)Zi\eta _3k_{xxxx}(0,0)Z`$ (76)
$`f_{xy}(x,y,t)\mathrm{\Theta }`$ $`=`$ $`i\eta _2k_{xxxx}(0,0)Z`$ (77)
$`f_{yy}(x,y,t)\mathrm{\Theta }`$ $`=`$ $`i\lambda k_{xx}(0,0)Zi\eta _1k_{xxxx}(0,0)Zi\eta _3k_{xxxx}(0,0)Z,`$ (78)
where $`k(xx^{},yy^{})=2D_0D(xx^{},yy^{}),k(0,0)=\frac{2D_0}{\pi \sigma ^2},k_{xx}(0,0)=k_{yy}(0,0)=\frac{4D_0}{\pi \sigma ^4}`$ and $`k_x(0,0)=k_y(0,0)=0`$. So we have

$`Z_t`$ $`=`$ $`i\gamma (t)\lambda Zi\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _1\mu _1}i\lambda {\displaystyle \frac{\alpha }{2}}Z_{\mu _2\mu _2}+i\alpha Z_{\mu _1x}+i\alpha Z_{\mu _2y}`$ (79)
$``$ $`i\alpha Z_{\eta _1}i\alpha Z_{\eta _3}+i\alpha \eta _1Z_{\eta _1\eta _1}+i\alpha \eta _1Z_{\eta _2\eta _2}`$ (80)
$`+`$ $`i\alpha \eta _2Z_{\eta _1\eta _2}+i\alpha \eta _2Z_{\eta _2\eta _3}+i\alpha \eta _3Z_{\eta _2\eta _2}+i\alpha \eta _3Z_{\eta _3\eta _3}`$ (81)
$``$ $`\lambda ^2k(0,0)Z+(\mu _1^2+\mu _2^22\lambda \eta _12\lambda \eta _3)k_{xx}(0,0)Z`$ (82)
$``$ $`(\eta _1^2+\eta _2^2+\eta _3^2+2\eta _1\eta _3)k_{xxxx}(0,0)Z.`$ (83)
Assuming statistical homogeneity $`(Z_x=0,Z_y=0)`$ and defining $`P(\stackrel{~}{h},u,v,w,s,q,t)`$ as the joint probability density function of $`\stackrel{~}{h},u,v,w,s`$ and $`q`$, one can construct the PDF as the Fourier transform of the generating function $`Z`$ with respect to $`\lambda ,\mu _1,\mu _2,\eta _1,\eta _2,\eta _3`$

$`P(\stackrel{~}{h},u,v,w,s,q,t)={\displaystyle }{\displaystyle \frac{d\lambda }{2\pi }}{\displaystyle \frac{d\mu _1}{2\pi }}{\displaystyle \frac{d\mu _2}{2\pi }}{\displaystyle \frac{d\eta _1}{2\pi }}{\displaystyle \frac{d\eta _2}{2\pi }}{\displaystyle \frac{d\eta _3}{2\pi }}\{`$ (84)
$`\mathrm{exp}(i\lambda \stackrel{~}{h}+i\mu _1u+i\mu _2v+i\eta _1w+i\eta _2s+i\eta _3q)Z(\lambda ,\mu _1,\mu _2,\eta _1,\eta _2,\eta _3,t)\}.`$ (85)
From the eqs. (43) and (44) the equation governing the evolution of $`P(\stackrel{~}{h},u,v,w,s,q,t)`$ can be derived, which is

$`P_t`$ $`=`$ $`\gamma (t)P_{\stackrel{~}{h}}+{\displaystyle \frac{\alpha }{2}}(u^2+v^2)P_{\stackrel{~}{h}}4\alpha wP4\alpha qP`$ (86)
$``$ $`\alpha w^2P_w\alpha q^2P_q\alpha s^2P_w\alpha s^2P_q\alpha wsP_s\alpha qsP_s`$ (87)
$`+`$ $`k(0,0)P_{\stackrel{~}{h}\stackrel{~}{h}}+2k_{xx}(0,0)P_{\stackrel{~}{h}w}+2k_{xx}(0,0)P_{\stackrel{~}{h}q}k_{xx}(0,0)P_{uu}k_{xx}(0,0)P_{vv}`$ (88)
$`+`$ $`k_{xxxx}(0,0)P_{ww}+k_{xxxx}(0,0)P_{ss}+k_{xxxx}(0,0)P_{qq}+2k_{xxxx}(0,0)P_{qw}.`$ (89)
From the eq.(45), it is easy to see that the moments $`\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_5}`$ satisfy the following equation

$`{\displaystyle \frac{d}{dt}}\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_5}=n_0\gamma (t)\stackrel{~}{h}^{n_01}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_5}{\displaystyle \frac{\alpha n_0}{2}}\stackrel{~}{h}^{n_01}u^{n_1+2}v^{n_2}w^{n_3}s^{n_4}q^{n_5}`$ (90)
$`{\displaystyle \frac{\alpha n_0}{2}}\stackrel{~}{h}^{n_01}u^{n_1}v^{n_2+2}w^{n_3}s^{n_4}q^{n_5}4\alpha \stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3+1}s^{n_4}q^{n_5}`$ (91)
$`4\alpha \stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_5+1}+\alpha (n_3+2)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3+1}s^{n_4}q^{n_5}`$ (92)
$`+\alpha n_3\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_31}s^{n_4+2}q^{n_5}+\alpha (n_4+1)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3+1}s^{n_4}q^{n_5}`$ (93)
$`+\alpha (n_4+1)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_5+1}+\alpha (n_5+2)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_5+1}`$ (94)
$`+n_5(n_51)k_{xxxx}(0,0)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_52}+2n_3n_5k_{xxxx}(0,0)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_31}s^{n_4}q^{n_51}`$ (95)
$`+2n_0n_3k_{xx}(0,0)\stackrel{~}{h}^{n_01}u^{n_1}v^{n_2}w^{n_31}s^{n_4}q^{n_5}+2n_0n_5k_{xx}(0,0)\stackrel{~}{h}^{n_01}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_51}`$ (96)
$`n_1(n_11)k_{xx}(0,0)\stackrel{~}{h}^{n_0}u^{n_12}v^{n_2}w^{n_3}s^{n_4}q^{n_5}n_2(n_21)k_{xx}(0,0)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_22}w^{n_3}s^{n_4}q^{n_5}`$ (97)
$`+n_3(n_31)k_{xxxx}(0,0)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_32}s^{n_4}q^{n_5}+n_4(n_41)k_{xxxx}(0,0)\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_42}q^{n_5}`$ (98)
$`+\alpha n_5\stackrel{~}{h}^{n_0}u^{n_1}v^{n_2}w^{n_3}s^{n_4+2}q^{n_51}+n_0(n_01)k(0,0)\stackrel{~}{h}^{n_02}u^{n_1}v^{n_2}w^{n_3}s^{n_4}q^{n_5}.`$ (99)
By substituting different sort of values for $`n_0,n_1,n_2,n_3,n_4,n_5`$ we can find some coupled differential equations for different moments. For example, if we take $`n_0=n_1=\mathrm{}=n_5=0`$ we have

$$<w>+<q>=0$$
(100)
or

$$<𝐮>=\mathrm{𝟎}$$
(101)
where this is the same as the statistical hemogeneity condition. For $`n_0=1,n_1=n_2=n_3=n_4=n_5=0`$ we find

$$<\stackrel{~}{h}w>+<\stackrel{~}{h}q>=<u^2><v^2>$$
(102)
and for $`n_0=0,n_1=1,n_2=n_3=n_4=n_5=0`$ we find

$$<uw>+<uq>=0,$$
(103)
and if we assume the statistical hemogeneity we have

$$<uw>=<uu_x>=\frac{1}{2}<u^2>_x=0.$$
(104)
So $`<uw>=0`$ and then $`<uq>=0`$, also for $`n_0=n_1=0,n_2=1,n_3=n_4=n_5=0`$ we find $`<vq>=0,<vw>=0`$.

For $`n_0=0,n_1=2,n_2=n_3=n_4=n_5=0`$ we find

$$\frac{d}{dt}<u^2>=\alpha <u^2w>\alpha <u^2q>2k_{xx}(0,0)$$
(105)
and also we have $`<u^2w>=\frac{1}{3}<u^3>_x=0`$ by statistical hemogeneity and also
$`<u^2q>=<u^2v>_y2<uvs>`$, $`<u^2v>_y=0`$, so $`<u^2q>=2<uvs>`$ so that we have

$$\frac{d}{dt}<u^2>=2\alpha <uvs>2k_{xx}(0,0).$$
(106)
The corresponding differential equation for $`<uvs>`$ is

$$\frac{d}{dt}<uvs>=0.$$
(107)
If we assume that at $`t=0`$ the surface is flat so all the moments at $`t=0`$ are zero,so that $`<uvs>=0`$ and then we find

$$<u^2>=2k_{xx}(0,0)t.$$
(108)
Similar calculations give

$$<v^2>=2k_{xx}(0,0)t.$$
(109)
By focusing on the differential equation of the PDF it is deduced that this equation is invariant under $`uu`$ and $`vv`$ which is a consequence of the inversion and reflection symmetry of the KPZ equation. So it will result in the vanishing of the $`u`$ and $`v`$ odd moments

$$u^{2k+1}=v^{2k+1}=0.$$
(110)
Using eq.(46) the even moments of $`u`$ and $`v`$ can be calculated. For example, for $`u^4`$ we have

$$\frac{d}{dt}u^4=\alpha u^4w\alpha u^4q12k_{xx}(0)u^2.$$
(111)
On the other hand $`u^4w`$ can be written as

$$u^4u_x=\frac{1}{5}u^5_x,$$
(112)
which is zero by homogeneity, also for $`u^4q`$ we have

$$u^4v_y=u^4v_y4u^3vs=4u^3vs$$
(113)
in which $`u^4v_y`$ is zero by homogeneity. The differential equation for $`u^3vs`$ is

$$\frac{d}{dt}u^3vs=0,$$
(114)
so $`u^3vs=0`$. Therefore, it is obtained

$$u^4=12k_{xx}^2(0,0)t^2,$$
(115)
and we have the same for $`v^4`$. Also it can be found that $`u^2v^2=4k_{xx}^2(0,0)t^2`$. By continuing the above method all the moments $`u^nv^m`$ can be deduced.

On the other hand, the above calculations show that all the mixed moments $`u^nv^ms`$ are zero. In the following some of the moments of $`u`$ and $`v`$ and their combinations have been given

$`u^6`$ $`=`$ $`v^6=120k_{xx}^3(0,0)t^3`$ (116)
$`u^2v^4`$ $`=`$ $`u^4v^2=24k_{xx}^3(0,0)t^3`$ (117)
$`u^8`$ $`=`$ $`v^8=1680k_{xx}^4(0,0)t^4`$ (118)
$`u^6v^2`$ $`=`$ $`u^2v^6=240k_{xx}^4(0,0)t^4`$ (119)
$`u^4v^4`$ $`=`$ $`144k_{xx}^4(0,0)t^4,etc.`$ (120)
We aim to calculate all the momentes of $`\stackrel{~}{h}=h\overline{h}`$. Putting $`n_0=2,n_1=n_2=\mathrm{}=0`$ in the eq.(46), we find

$$\frac{d}{dt}\stackrel{~}{h}^2=\alpha \stackrel{~}{h}u^2\alpha \stackrel{~}{h}v^2\alpha \stackrel{~}{h}^2w\alpha \stackrel{~}{h}^2q+2k(0,0).$$
(121)
We know that $`\stackrel{~}{h}^2w=\stackrel{~}{h}^2u_x2\stackrel{~}{h}u^2`$ and because of the homogeneity $`\stackrel{~}{h}^2u_x=0`$, so we have

$$\stackrel{~}{h}^2w=2\stackrel{~}{h}u^2,$$
(122)
and also

$$\stackrel{~}{h}^2q=2\stackrel{~}{h}v^2,$$
(123)
then

$$\frac{d}{dt}\stackrel{~}{h}^2=\alpha (\stackrel{~}{h}u^2+\stackrel{~}{h}v^2)+2k(0,0).$$
(124)
As it can be seen from eq.(71) we need the moments $`\stackrel{~}{h}u^2`$ and $`\stackrel{~}{h}v^2`$ to find $`\stackrel{~}{h}^2`$ , which the related differential equations are

$`{\displaystyle \frac{d}{dt}}\stackrel{~}{h}u^2`$ $`=`$ $`\gamma (t)u^2{\displaystyle \frac{\alpha }{2}}u^4{\displaystyle \frac{\alpha }{2}}u^2v^2\alpha \stackrel{~}{h}u^2w\alpha \stackrel{~}{h}u^2q`$ (125)
$`{\displaystyle \frac{d}{dt}}\stackrel{~}{h}v^2`$ $`=`$ $`\gamma (t)v^2{\displaystyle \frac{\alpha }{2}}v^4{\displaystyle \frac{\alpha }{2}}u^2v^2\alpha \stackrel{~}{h}v^2w\alpha \stackrel{~}{h}v^2q.`$ (126)
By using of the statistical homogeneity the last two terms of the above equations can be converted as

$`\stackrel{~}{h}u^2w`$ $`=`$ $`{\displaystyle \frac{1}{3}}u^4`$ (127)
$`\stackrel{~}{h}v^2q`$ $`=`$ $`{\displaystyle \frac{1}{3}}v^4`$ (128)
$`\stackrel{~}{h}u^2q`$ $`=`$ $`u^2v^22\stackrel{~}{h}uvs`$ (129)
$`\stackrel{~}{h}v^2w`$ $`=`$ $`u^2v^22\stackrel{~}{h}uvs.`$ (130)
The above relations result in

$`{\displaystyle \frac{d}{dt}}(\stackrel{~}{h}u^2+\stackrel{~}{h}v^2`$ $`=`$ $`\gamma (t)(u^2+v^2){\displaystyle \frac{\alpha }{6}}(u^4+v^4)`$ (131)
$`+`$ $`\alpha u^2v^2+4\alpha \stackrel{~}{h}uvs.`$ (132)
Obtaining the differential equation for $`\stackrel{~}{h}uvs`$ we get

$$\frac{d}{dt}\stackrel{~}{h}uvs=0,$$
(133)
which results in $`\stackrel{~}{h}uvs=0`$. Also it is easy to see that all the moments $`\stackrel{~}{h}^nu^mv^ps`$ are zero. This fantastic result helps us to find all the $`\stackrel{~}{h}^n`$ moments. Now by substituting $`\gamma (t)=\overline{h}_t=\frac{\alpha }{2}(u^2+v^2)=2\alpha k_{xx}(0,0)t`$ , $`u^4`$ , $`v^4`$ and $`u^2v^2`$ in eq.(75) we obtain

$$\stackrel{~}{h}u^2+\stackrel{~}{h}v^2=\frac{8}{3}\alpha k_{xx}^2(0,0)t^3,$$
(134)
which finally gives

$$\stackrel{~}{h}^2=\frac{2}{3}\alpha ^2k_{xx}^2(0,0)t^4+2k(0,0)t.$$
(135)
Now we begin to calculate the moment $`\stackrel{~}{h}^3`$. Inserting $`n_0=3,n_1=n_2=\mathrm{}=n_5=0`$ in eq.(46) we get

$$\frac{d}{dt}\stackrel{~}{h}^3=3\gamma (t)\stackrel{~}{h}^2\frac{3}{2}\alpha (\stackrel{~}{h}^2u^2+\stackrel{~}{h}^2v^2)\alpha \stackrel{~}{h}^3w\alpha \stackrel{~}{h}^3q,$$
(136)
and again by using statisiacal homogeneity it can be shown that

$`\stackrel{~}{h}^3w`$ $`=`$ $`3\stackrel{~}{h}^2u^2`$ (137)
$`\stackrel{~}{h}^3q`$ $`=`$ $`3\stackrel{~}{h}^2v^2.`$ (138)
So we have

$$\frac{d}{dt}\stackrel{~}{h}^3=3\gamma (t)\stackrel{~}{h}^2+\frac{3}{2}\alpha (\stackrel{~}{h}^2u^2+\stackrel{~}{h}^2v^2).$$
(139)
To calculate $`\stackrel{~}{h}^3`$ we need $`\stackrel{~}{h}^2`$ , $`(\stackrel{~}{h}^2u^2+\stackrel{~}{h}^2v^2)`$. $`\stackrel{~}{h}^2`$ has been calculated above, so we will obtain $`\stackrel{~}{h}^2u^2+\stackrel{~}{h}^2v^2`$ using the corresponding differential equation as follows,

$`{\displaystyle \frac{d}{dt}}(\stackrel{~}{h}^2u^2+\stackrel{~}{h}^2v^2`$ $`=`$ $`2\gamma (t)(\stackrel{~}{h}u^2+\stackrel{~}{h}v^2)\alpha (\stackrel{~}{h}u^4+\stackrel{~}{h}v^4)2\alpha \stackrel{~}{h}u^2v^2\alpha (\stackrel{~}{h}^2u^2w+\stackrel{~}{h}^2v^2w)`$ (140)
$``$ $`\alpha (\stackrel{~}{h}^2u^2q+\stackrel{~}{h}^2v^2q)+2k(0,0)(u^2+v^2)4k_{xx}(0,0)\stackrel{~}{h}^2.`$ (141)
As before, we easily get

$`\stackrel{~}{h}^2u^2w`$ $`=`$ $`{\displaystyle \frac{2}{3}}\stackrel{~}{h}u^4`$ (142)
$`\stackrel{~}{h}^2v^2q`$ $`=`$ $`{\displaystyle \frac{2}{3}}\stackrel{~}{h}v^4`$ (143)
$`\stackrel{~}{h}^2u^2q`$ $`=`$ $`2\stackrel{~}{h}u^2v^22\stackrel{~}{h}^2uvs`$ (144)
$`\stackrel{~}{h}^2v^2w`$ $`=`$ $`2\stackrel{~}{h}u^2v^22\stackrel{~}{h}^2uvs.`$ (145)
As discussed before, it is easy to show that the moment $`\stackrel{~}{h}^2uvs`$ is zero. To prove this we write the corresponding differential equation

$$\frac{d}{dt}\stackrel{~}{h}^2uvs=\alpha (\stackrel{~}{h}u^3vs+\stackrel{~}{h}uv^3s),$$
(146)
and again by trying to write the differential eqautions for $`\stackrel{~}{h}u^3vs`$ and $`\stackrel{~}{h}uv^3s`$ we obtain

$`{\displaystyle \frac{d}{dt}}\stackrel{~}{h}u^3vs`$ $`=`$ $`0`$ (147)
$`{\displaystyle \frac{d}{dt}}\stackrel{~}{h}uv^3s`$ $`=`$ $`0`$ (148)
which results in $`\stackrel{~}{h}u^3vs=\stackrel{~}{h}uv^3s=0`$, therefore $`\stackrel{~}{h}^2uvs=0`$. Now $`\stackrel{~}{h}u^4`$ , $`\stackrel{~}{h}v^4`$ and $`\stackrel{~}{h}u^2v^2`$ should be found. The relating differential equation for $`\stackrel{~}{h}u^4+\stackrel{~}{h}v^4`$ is

$`{\displaystyle \frac{d}{dt}}(\stackrel{~}{h}u^4+\stackrel{~}{h}v^4`$ $`=`$ $`\gamma (t)(u^4+v^4){\displaystyle \frac{\alpha }{2}}(u^6+v^6){\displaystyle \frac{\alpha }{2}}(u^4v^2+u^2v^4)`$ (149)
$``$ $`\alpha (\stackrel{~}{h}u^4w+\stackrel{~}{h}v^4w)\alpha (\stackrel{~}{h}u^4q+\stackrel{~}{h}v^4q)12k_{xx}(0,0)(\stackrel{~}{h}u^2+\stackrel{~}{h}v^2).`$ (150)
As before the following identities are held

$`\stackrel{~}{h}u^4w`$ $`=`$ $`{\displaystyle \frac{1}{5}}u^6`$ (151)
$`\stackrel{~}{h}u^4q`$ $`=`$ $`u^4v^2`$ (152)
$`\stackrel{~}{h}v^4q`$ $`=`$ $`{\displaystyle \frac{1}{5}}v^6`$ (153)
$`\stackrel{~}{h}v^4w`$ $`=`$ $`u^2v^4,`$ (154)
so we have

$`{\displaystyle \frac{d}{dt}}(\stackrel{~}{h}u^4+\stackrel{~}{h}v^4`$ $`=`$ $`\gamma (t)(u^4+v^4){\displaystyle \frac{3\alpha }{10}}(u^6+v^6)+{\displaystyle \frac{\alpha }{2}}(u^4v^2+u^2v^4)`$ (155)
$``$ $`12k_{xx}(0,0)(\stackrel{~}{h}u^2+\stackrel{~}{h}v^2).`$ (156)
Substituting the expressions for $`\gamma (t)`$ ,$`u^4`$ , $`v^4`$ , $`u^6`$ , $`v^6`$ , $`u^4v^2`$, $`u^2v^4`$ and $`(\stackrel{~}{h}u^2+\stackrel{~}{h}v^2)`$

we find

$$\stackrel{~}{h}u^4+\stackrel{~}{h}v^4=32\alpha k_{xx}(0,0)^3t^4.$$
(157)
The corresponding differential equation for $`\stackrel{~}{h}u^2v^2`$ is

$`{\displaystyle \frac{d}{dt}}(\stackrel{~}{h}u^2v^2`$ $`=`$ $`\gamma (t)(u^2v^2){\displaystyle \frac{\alpha }{2}}(u^4v^2+u^2v^4)`$ (158)
$``$ $`\alpha (\stackrel{~}{h}u^2v^2w+\stackrel{~}{h}u^2v^2q)2k_{xx}(0,0)(\stackrel{~}{h}u^2+\stackrel{~}{h}v^2).`$ (159)
As before the following identities are held

$`\stackrel{~}{h}u^2v^2w`$ $`=`$ $`{\displaystyle \frac{1}{3}}u^4v^2`$ (160)
$`\stackrel{~}{h}u^2v^2q`$ $`=`$ $`{\displaystyle \frac{1}{3}}u^2v^4,`$ (161)
then by inserting the known moments we find

$$\stackrel{~}{h}u^2v^2=\frac{16}{3}\alpha k_{xx}^3(0,0)t^4.$$
(162)
The above results get

$$\stackrel{~}{h}^2u^2+\stackrel{~}{h}^2v^2=\frac{8}{3}\alpha ^2k_{xx}^3(0,0)t^58k(0,0)k_{xx}(0,0)t^2,$$
(163)
which finally results in

$$\stackrel{~}{h}^3=\frac{48}{45}\alpha ^3k_{xx}^3(0,0)t^6.$$
(164)
By continuing the above proccedure all the $`\stackrel{~}{h}^n`$ moments can be derived. Some of these moments are listed bellow

$`\stackrel{~}{h}^4`$ $`=`$ $`{\displaystyle \frac{44}{35}}\alpha ^4k_{xx}^4(0,0)t^88\alpha ^2k(0,0)k_{xx}^2(0,0)t^5+12k^2(0,0)t^2`$ (165)
$`\stackrel{~}{h}^5`$ $`=`$ $`{\displaystyle \frac{1216}{945}}\alpha ^5k_{xx}^5(0,0)t^{10}{\displaystyle \frac{64}{3}}\alpha ^3k(0,0)k_{xx}^3(0,0)t^7,`$ (166)
which are the same as the results that we derived by expanding the generating function. Results of this appendix show that any moment containing the first power of $`s`$ vanishes. Indeed, one can prove the following identity

$$se^{i(\lambda \stackrel{~}{h}+\mu _1u+\mu _2v)}=0.$$
(167)
By expanding the exponential in the above expression one finds

$$se^{i(\lambda \stackrel{~}{h}+\mu _1u+\mu _2v)}=\underset{n,m,p}{}\frac{i^{(n+m+p)}\lambda ^n\mu _1^m\mu _2^p}{n!m!p!}s\stackrel{~}{h}^nu^mv^p.$$
(168)
Now setting $`n_3=n_5=0`$ and $`n_4=1`$ in eq.(46) we get the following equation for $`s\stackrel{~}{h}^nu^mv^p`$

$`{\displaystyle \frac{d}{dt}}\stackrel{~}{h}^nu^mv^ps=`$ $``$ $`n\gamma (t)\stackrel{~}{h}^{n1}u^mv^ps{\displaystyle \frac{\alpha n}{2}}\stackrel{~}{h}^{n1}u^{m+2}v^ps{\displaystyle \frac{\alpha n}{2}}\stackrel{~}{h}^{n1}u^mv^{p+2}s`$ (169)
$`+`$ $`n(n1)k_{xx}(0,0)\stackrel{~}{h}^{n2}u^mv^psm(m1)k_{xx}(0,0)\stackrel{~}{h}^nu^{m2}v^ps`$ (171)
$``$ $`p(p1)k_{xx}(0,0)\stackrel{~}{h}^nu^mv^{p2}s.`$ (173)
Starting from eqs.(54),(61) and (76) all the $`s\stackrel{~}{h}^nu^mv^p`$ moments can be evaluated and it can be shown that, with flat initial condition, they are equal to zero. So we conclude $`s\mathrm{\Theta }=u_y\mathrm{\Theta }=v_x\mathrm{\Theta }=0`$.


text:
# References

Supersymmetric generalizations of integrable systems have been a field of considerable interest since the KdV equation was incorporated in a supersymmetric framework . Apart from being a nonlinear equation of widespread applicability, it has added significance as equations of the KdV type arise in the context of 2D quantum gravity . Moreover its supersymmetric extensions are important, since it is believed that a fundamental physical theory must possess supersymmetry. Various aspects of integrability of supersymmetric equations have been extensively studied mainly by extending the Lax formalism to the supersymmetric situation. The biHamiltonian structures of the superKdV as well as the superKP hierarchies have been obtained by suitably formulating the Gelfand Dikii method for the supersymmetric systems . Many of the techniques used in the bosonic integrable systems have been extended to the supersymmetric version of the theory. Mention may be made of the Backlund tranformations , the Painleve analysis , the Darboux transformations  and the $`\tau `$-functions  vis a vis the soliton solutions which have been generalized to the supersymmetric framework. However, a powerful and general method of generating soliton solutions in bosonic nonlinear differential equations, namely the Hirota bilinear method  has been applied to supersymmetric extensions of only a few nonlinear systems; these have been mainly $`N=1`$ supersymmetric equations . But the bilinear method for $`N=2`$ supersymmetric equations are yet to be studied.

A supersymmetric system of special interest is the $`N=2`$ KdV equation formulated by Inami and Kanno  from the non-standard superLax representation. It was later shown that this system could be obtained from the reduction of a $`N=2`$ superKP hierarchy  arising out of the non-standard flow equation. Even though integrability of this system has been established from the existence of biHamiltonian structures , its bilinear form as well as its soliton solutions are known only when certain conditions are imposed on the superfields . In this letter, an attempt is made to bilinearize $`N=2`$ supersymmetric systems, embedded in $`N=1`$ superspace, following Hirota’s method and to obtain soliton solutions. This will be a generalization of the earlier work  in the sense that no constraints on the superfields have been imposed. As an example, the supersymmetric generalization of Hirota bilinear method is considered for $`N=2`$ supersymmetric KdV equation of Inami and Kanno. Specifically, soliton solutions of the Inami Kanno equation are obtained by casting it in the bilinear form and existence of $`N`$ soliton solutions are shown to exist thereby establishing the integrability of the system.

The Inami Kanno equation  is given by

$$_tu_1+u_1^{[6]}+3\left(u_1^{[1]}u_0\right)^{[2]}\frac{1}{2}\left(u_1^3\right)^{[2]}=0$$
(1)
$$_tu_0+u_0^{[6]}3\left(u_0u_0^{[1]}\right)^{[2]}\frac{3}{2}\left(u_0u_1^2\right)^{[2]}+3\left(u_0u_1^{[2]}\right)^{[2]}=0$$
(2)
where $`[n]`$ denotes the $`n`$th derivative with respect to the $`D`$ operator, defined by

$$D=\frac{}{\theta }+\theta \frac{}{x}$$
(3)
and $`u_1`$ and $`u_0`$ are superfields of conformal spin $`1`$ and $`3/2`$ respectively. It is clear from the form of the equations that they are two coupled $`N=1`$ supersymmetric equations which in the limits $`u_1=0`$ and $`u_0=0`$ reduce to the $`N=1`$ KdV equation and the $`N=1`$ mKdV equation  respectively. Interestingly, in yet another limit, namely $`u_0=Du_1`$ both the equations (1,2) acquire an identical form and ultimately reduce to the $`N=1`$ mKdV equation. All these reductions are bilinearizable and $`N`$ soliton solutions exist for each of them . On the other hand, the bosonic limit of (1,2) gives rise to a set of coupled equations of the form

$$_tu_1^b+^3u_1^b\frac{1}{2}u_1^{b^3}=0$$
(4)
and

$$_tu_0^b+^3u_0^b3u_0^{b^2}\frac{3}{2}\left(u_0^bu_1^{b^2}\right)+3\left[u_0^b(u_1^b)\right]=0$$
(5)
where the super fields are written in the component forms as

$$u_1=u_1^b+\theta u_1^f$$
(6)
and

$$u_0=u_0^f+\theta u_0^b$$
(7)
The first equation, (4) is nothing but the mKdV equation, whereas the second one, (5) is a coupled one and reduces to the KdV equation if $`u_1^b`$ is set to zero. The bosonic equations can be shown to be integrable and possess $`N`$ soliton solutions.

In order to cast the supersymmetric equations (1,2) in the bilinear form, we transform the superfields as

$$u_1=2D^2\mathrm{log}\frac{\tau _1}{\tau _2}$$
(8)
for $`u_1`$ and

$$u_0=2D^3\mathrm{log}\tau _1$$
(9)
for $`u_0`$. $`\tau _1`$ and $`\tau _2`$ are independent bosonic superfields. Notice that in (9), $`u_0`$ depends only the superfield $`\tau _1`$, whereas the other superfield $`u_1`$ is a function of both $`\tau _1`$ and $`\tau _2`$.

In terms of the supersymmetric Hirota derivative $`𝐒`$ , defined by

$$\mathrm{𝐒𝐃}_x^nf.g=(D_{\theta _1}D_{\theta _2})(_{x_1}_{x_2})^nf(x_1,\theta _1)g(x_2,\theta _2)|_{\begin{array}{c}x_1=x_2=x\\ \theta _1=\theta _2=\theta \end{array}}$$
(10)
where $`𝐃`$ is ordinary (or bosonic) Hirota derivative, the $`N=2`$ super KdV equations can be cast in the following bilinear form:

$$(\mathrm{𝐒𝐃}_t+\mathrm{𝐒𝐃}_x^3)(\tau _1.\tau _1)=0$$
(11)
$$(\mathrm{𝐒𝐃}_t+\mathrm{𝐒𝐃}_x^3)(\tau _2.\tau _2)=0$$
(12)
$$𝐃_x^2(\tau _1.\tau _2)=0$$
(13)
The bilinear forms in (11,12,13), indeed, ensure the existence of one and two soliton solutions of $`N=2`$ KdV equations (1,2).

For one soliton solution, we use the the following forms of the $`\tau `$-functions:

$$\tau _1=1+\alpha e^\eta $$
(14)
and

$$\tau _2=1+\beta e^\eta $$
(15)
$`\alpha `$ and $`\beta `$ in (14,15) are nonzero bosonic parameters. The parameters $`\alpha `$ and $`\beta `$ can be determined from the initial value conditions. In (14,15),

$$\eta =kx+\omega t+\zeta \theta $$
(16)
where $`k`$ and $`\omega `$ are the bosonic parameters and $`\zeta `$ is the Grassmann odd parameter.

Substituting (14) and (15) back into third equation of the Hirota form of the superKdV equation (13), we find the non-trivial solutions exist provided

$$\beta =\alpha .$$
(17)
The dispersion relation, however, follows from (11,12) as

$$\omega +k^3=0$$
(18)
which is identical to that for the bosonic KdV or mKdV equations. In obtaining the dispersion relation the following property of the superHirota operator has been used.

$$\mathrm{𝐒𝐃}^n(e^{\eta _1}.e^{\eta _2})=(k_1k_2)^n[(\zeta _1\zeta _2)+\theta (k_1k_2)]e^{\eta _1+\eta _2}$$
(19)
The fermionic parameter, $`\zeta `$, however, remains arbitrary at the one soliton level.

The explicit forms of the the soliton solutions for the superfields $`u_1`$ and $`u_0`$ may be found by substituting the $`\tau `$\- functions in (8) and (9) respectively. Explicitly

$$u_1=(2k)\text{cosech}(\varphi +\gamma _0)\theta (2k\zeta )\text{cosh}(\varphi +\gamma _0)\text{cosech}^2(\varphi +\gamma _0)$$
(20)
and

$$u_0=\frac{1}{2}k\zeta \text{sech}^2\frac{1}{2}(\varphi +\gamma _0)+\theta \frac{k^2}{2}\text{sech}^2\frac{1}{2}(\varphi +\gamma _0)$$
(21)
where we have chosen $`\varphi =kxk^3t`$ and $`\alpha =\beta =e^{\gamma _0}`$, $`\gamma _0`$ being nonzero, real parameter. The bosonic and fermionic components of the superfields are easily found by comparing (20) and (21) with (6) and (7). It is observed that the soliton solution of $`u_1^b`$ becomes isomorphic to the mKdV soliton and $`u_0^b`$ to the KdV soliton as expected.

If we choose $`\alpha `$ to be negative in (14), the solutions for $`u_1`$ and $`u_0`$ become

$$u_1=(2k)\text{cosech}(\varphi +\gamma _0)+\theta (2k\zeta )\text{cosh}(\varphi +\gamma _0)\text{cosech}^2(\varphi +\gamma _0)$$
(22)
and

$$u_0=\frac{1}{2}k\zeta \text{cosech}^2\frac{1}{2}(\varphi +\gamma _0)\theta \frac{k^2}{2}\text{cosech}^2\frac{1}{2}(\varphi +\gamma _0)$$
(23)
Interestingly, the solutions (20,21) and (22,23) are connected through the transformations

$$u_1=u_1$$
(24)
and

$$u_0=u_0Du_1$$
(25)
On the other hand, at the equation level, the Inami Kanno equations remain invariant with respect to the transformations (24,25). Evidently, $`\alpha \alpha `$ is a symmetry of the the equations of motion. At the $`\tau `$ functions level (8,9), in fact, this symmetry is manifested in the tranformation $`\tau _1\tau _2`$ and vice versa and thus will be observed at the $`N`$ soliton solutions.

The existence of one soliton solution, however, does not ensure the exact integrability in the Hirota formalism. In order that integrability be established, at least three soliton solution must exist. The one soliton solution fixes the arbitrary parameters in the trial solutions. On the other hand, the interaction term is determined by the two soliton solutions, and its form remains unchanged for all higher solitons if the system is truly integrable. In fact, if the equations of motion can be cast in the bilinear form, one and two soliton solutions can be always determined. It is thus necessary to show the existence of at least three soliton solution, which checks the consistency of the solution.

The $`\tau _1`$ for the $`N`$ soliton of the $`N=2`$ superKdV may be written as

$$\tau _1=\underset{\mu _i=0,1}{}exp\left(\underset{i,j=1}{\overset{N}{}}\varphi (i,j)\mu _i\mu _j+\underset{i=1}{\overset{N}{}}\mu _i(\eta _i+\mathrm{log}\alpha _i)\right)(i<j)$$
(26)
where $`A_{ij}=e^{\varphi (i,j)}`$. For $`\tau _2`$ we replace $`\alpha _i`$ by $`\beta _i`$ and $`A_{ij}`$ by $`B_{ij}`$. In particular, for the two soliton solution of the superKdV equations $`\tau `$ functions may be written in the form

$$\tau _1=1+\alpha _1e^{\eta _1}+\alpha _2e^{\eta _2}+\alpha _1\alpha _2A_{12}e^{\eta _1+\eta _2}$$
(27)
and

$$\tau _2=1+\beta _1e^{\eta _1}+\beta _2e^{\eta _2}+\beta _1\beta _2B_{12}e^{\eta _1+\eta _2}$$
(28)
where

$$\eta _1=k_1x+\omega _1t+\zeta _1\theta $$
(29)
and

$$\eta _2=k_2x+\omega _2t+\zeta _2\theta $$
(30)
In (29,30) the parameters $`k_1`$, $`k_2`$, $`\omega _1`$ and $`\omega _2`$ are bosonic, while $`\zeta _1`$ and $`\zeta _2`$ are fermionic ones, as before. If we consider $`\alpha _1=\beta _1`$, $`\alpha _2=\beta _2`$,

$$A_{12}=B_{12}=\frac{(k_1k_2)^2}{(k_1+k_2)^2}$$
(31)
and

$$k_1\zeta _2=k_2\zeta _1$$
(32)
the two soliton solutions in (27,28) satisfy the the bilinear forms (11,12,13) leading to the dispersion relations

$$\omega _1+k_1^3=0$$
(33)
and

$$\omega _2+k_2^3=0$$
(34)
The condition (32) relating fermionic and bosonic parameters is necessary to obtain $`A_{12}`$ and $`B_{12}`$, given in (31) consistently. It will also be seen that the constraint on the parameters is crucial in demonstrating the existence of three soliton solutions.

The explicit forms of $`\tau _1`$ and $`\tau _2`$ for the three soliton solution are

$`\tau _1=1+\alpha _1e^{\eta _1}+\alpha _2e^{\eta _2}+\alpha _3e^{\eta _3}+\alpha _1\alpha _2A_{12}e^{\eta _1+\eta _2}+\alpha _1\alpha _3A_{13}e^{\eta _1+\eta _3}`$
$`+\alpha _2\alpha _3A_{23}e^{\eta _2+\eta _3}+\alpha _1\alpha _2\alpha _3A_{12}A_{13}A_{23}e^{\eta _1+\eta _2+\eta _3}`$ (35)
and

$`\tau _2=1+\beta _1e^{\eta _1}+\beta _2e^{\eta _2}+\beta _3e^{\eta _3}+\beta _1\beta _2B_{12}e^{\eta _1+\eta _2}+\beta _1\beta _3B_{13}e^{\eta _1+\eta _3}`$
$`+\beta _2\beta _3B_{23}e^{\eta _2+\eta _3}+\beta _1\beta _2\beta _3B_{12}B_{13}B_{23}e^{\eta _1+\eta _2+\eta _3}`$ (36)
where

$$\eta _i=k_{ix}x+\omega _it+\zeta _i\theta (i=1,2,3)$$
(37)
Notice that the three soliton solution does not contain any new unknown parameter; it is expressed in terms of the parameters of the two soliton solutions only. These forms of three soliton solution give rise to a set of nontrivial relations among the parameters, which determine the consistency of the solutions. In the three soliton solution we find the interaction term from (11) or (12), in the same way as for two solitons. Explicitly,

$$A_{ij}=B_{ij}=\frac{(k_ik_j)^2}{(k_i+k_j)^2}(i,j=1,2,3;ij)$$
(38)
This, in turn, imposes the the following constraint on the fermionic parameters:

$$k_i\zeta _j=k_j\zeta _i(i,j=1,2,3;ij)$$
(39)
The interaction term $`A_{ij}`$ can also be found from the third bilinear form (13) and is identical to (38). It is found by direct substitution that the three soliton solutions (35) and (36) are consistent with conditions derived from the one and two soliton solutions. This proves that the $`N=2`$ super KdV equation of Inami and Kanno is also integrable from the Hirota point of view.

In conclusion, this is a first attempt to bilinearize $`N=2`$ supersymmetric equation and to obtain $`N`$ soliton solutions. The example considered here is the $`N=2`$ KdV equation of Inami and Kanno which has two interesting $`N=1`$ limits. If $`u_1`$ is set to zero, it reduces to $`N=1`$ KdV equation, whereas by choosing the limit $`u_0=0`$ or $`u_0=Du_1`$, it becomes $`N=1`$ mKdV equation. The fermionic parameter $`\zeta `$ plays a trivial role at the one soliton solution level, but at the higher soliton level it is intimately related with the bosonic parameters and eventually the conditions become essential to show the existence of $`N`$ soliton solutions vis a vis the integrability of the system.

SG would like to thank DST, Govt. of India for financial support under the project no. 100/(IFD)/2066/2000-2002.