text:
# Non-Gaussian Statistics of the Soliton Timing Jitter due to Amplifier Noise

## Abstract

Based on the first-order perturbation theory of soliton, the Gordon-Haus timing jitter induced by amplifier noise is found to be non-Gaussian distributed. Compared with Gaussian distribution given by the linearized perturbation theory, both frequency and timing jitter have larger tail probability. The timing jitter has a larger discrepancy to Gaussian distribution than that of frequency jitter.

fiber soliton, timing jitter, Gordon-Haus effect

The Gordon-Haus timing jitter in fiber soliton due to amplifier noise is usually assumed to be Gaussian distributed Gordon and Haus (1986); Moore et al. (2003); Iannone et al. (1998) when the first-order perturbation theory of soliton Kivshar and Malomed (1989); Kaup (1990); Georges (1995a) is used. Previous works showed that the non-Gaussian timing jitter is induced by soliton interactions Menyuk (1995); Georges (1995b, 1996) and regeneration Leclerc et al. (1999); Falkovich et al. (2001) but not by amplifier noise alone. When the first-order perturbation of soliton is linearized Gordon and Haus (1986); Iannone et al. (1998), the Gordon-Haus timing jitter is indeed Gaussian distributed. However, if the ordinary differential equations from the first-order perturbation is not linearized, as shown later, the amplitude, frequency and timing jitters are all non-Gaussian distributed.

In this letter, the amplitude jitter is found to be noncentral chi-square ($`\chi ^2`$) distributed, confirming the previous simulation results Moore et al. (2003) and the non-Gaussian statistics Falkovich et al. (2001). Both the frequency and timing jitters are also non-Gaussian distributed.

From the first-order perturbation theory of soliton Kivshar and Malomed (1989); Kaup (1990); Georges (1995a), with amplifier noise alone, the soliton parameters are evolved according to the following ordinary (or stochastic) differential equations:

$`{\displaystyle \frac{dA}{d\zeta }}`$ $`=`$ $`n_A(\zeta ),`$ (1)
$`{\displaystyle \frac{d\mathrm{\Omega }}{d\zeta }}`$ $`=`$ $`n_\mathrm{\Omega }(\zeta ),`$ (2)
$`{\displaystyle \frac{dT}{d\zeta }}`$ $`=`$ $`\mathrm{\Omega }+n_T(\zeta ).`$ (3)
Not used here, the phase perturbation is not shown in (1)-(3). All noise terms of $`n_A(\zeta ),n_\mathrm{\Omega }(\zeta )`$, and $`n_T(\zeta )`$ are independent Gaussian processes with autocorrelation of Georges (1995a, 1996)

$`E\{n_A(\zeta _1)n_A(\zeta _2)\}`$ $`=`$ $`A\sigma _n^2\delta (\zeta _1\zeta _2),`$ (4)
$`E\{n_\mathrm{\Omega }(\zeta _1)n_\mathrm{\Omega }(\zeta _2)\}`$ $`=`$ $`{\displaystyle \frac{A}{3}}\sigma _n^2\delta (\zeta _1\zeta _2),`$ (5)
$`E\{n_T(\zeta _1)n_T(\zeta _2)\}`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{12A^3}}\sigma _n^2\delta (\zeta _1\zeta _2),`$ (6)
where $`\sigma _n^2`$ is the noise variance and $`E\{\}`$ denotes expectation. From (4) to (6), the variances depend on the amplitude of $`A(\zeta )`$ and amplitude jitter transfers to both frequency and timing jitter. If the amplitude in the variances of (4) to (6) is assumed to be constant \[$`A(\zeta )A(0)=A_0`$\] as a first-order approximation, amplitude, frequency and timing jitters are indeed Gaussian distributed.

The inclusion of amplitude jitter in (4) to (6) is still within the equations of the first-order soliton perturbation theory Kivshar and Malomed (1989); Kaup (1990); Georges (1995a). The nonlinear first-order perturbation can be interpreted the repeated usage of the linearized first-order perturbation Iannone et al. (1998).

Based on (1) and (4), the stochastic differential equation (SDE) of the amplitude jitter is $`dA=\sqrt{A}dw_A`$ where $`w_A`$ is a Wiener process with autocorrelation of $`E\{w_A(\zeta _1)w_A(\zeta _2)\}=\sigma _n^2\mathrm{min}(\zeta _1,\zeta _2)`$. With an initial value of $`A(0)=A_0`$, the amplitude jitter is

$$A(\zeta )=\left(A_0^{1/2}+\frac{w_A}{2}\right)^2$$
(7)
as a noncentral $`\chi ^2`$-distributed random process Humblet and Azizog̃lu (1991); Shum and Ghafouri-Shiraz (1996); Moore et al. (2003) with variance parameter of $`\frac{1}{4}\sigma _n^2\zeta `$. The amplitude jitter of (7) is the solution of a Stranovich but not Ito SDE Gardiner (1985).

Using (2) and (5), the frequency jitter is

$$\mathrm{\Omega }(\zeta )=_0^\zeta \left(A_0^{1/2}+\frac{w_A(\zeta _1)}{2}\right)𝑑w_\mathrm{\Omega },$$
(8)
where $`w_\mathrm{\Omega }`$ is a Wiener process with autocorrelation of $`E\{w_\mathrm{\Omega }(\zeta _1)w_\mathrm{\Omega }(\zeta _2)\}=\frac{1}{3}\sigma _n^2\mathrm{min}(\zeta _1,\zeta _2)`$ and independent of the Wiener process of $`w_A`$.

The timing jitter equation of (3) has two terms, the first-term of $`\mathrm{\Omega }`$ gives the Gordon-Haus timing jitter increasing with $`\zeta ^3`$ and second term of $`n_T(\zeta )`$ is just the projection of amplifier noise into timing jitter Gordon and Haus (1986). The first term of Gordon-Haus timing jitter is far more interesting than the second term. The SDE of $`dT_{\mathrm{GH}}=\mathrm{\Omega }d\zeta `$ has a solution of

$$T_{\mathrm{GH}}(\zeta )=_0^\zeta (\zeta \zeta _1)\left(A_0^{1/2}+\frac{w_A(\zeta _1)}{2}\right)𝑑w_\mathrm{\Omega }.$$
(9)
Similar to option pricing with stochastic volatility Stein and Stein (1991), the characteristic functions of the frequency $`\mathrm{\Omega }(\zeta )`$ and Gordon-Haus timing jitter $`T_{\mathrm{GH}}(\zeta )`$ are

$`\mathrm{\Psi }_{\mathrm{\Omega }(\zeta )}(\nu )=G_1\left({\displaystyle \frac{\nu ^2\sigma _n^2}{6}}\right),`$ (10)
$`\mathrm{\Psi }_{T_{\mathrm{GH}}(\zeta )}(\nu )=G_2\left({\displaystyle \frac{\nu ^2\sigma _n^2}{6}}\right),`$ (11)
and

$`G_1(\lambda )`$ $`=`$ $`E\left\{e^{\lambda _0^\zeta \left(A_0^{1/2}+\frac{w_A(\zeta _1)}{2}\right)^2𝑑\zeta _1}\right\},`$ (12)
$`G_2(\lambda )`$ $`=`$ $`E\left\{e^{\lambda _0^\zeta (\zeta \zeta _1)^2\left(A_0^{1/2}+\frac{w_A(\zeta _1)}{2}\right)^2𝑑\zeta _1}\right\}`$ (13)
where $`G_1(\lambda )`$ and $`G_2(\lambda )`$ are the moment generating functions of $`_0^\zeta \left(A_0^{1/2}+\frac{1}{2}w_A\right)^2𝑑\zeta _1`$ and $`_0^\zeta (\zeta \zeta _1)^2\left(A_0^{1/2}+\frac{1}{2}w_A\right)^2𝑑\zeta _1`$, respectively.

Based on the Cameron-Martin integral Cameron and Martin (1945), we get

$$G_1(\lambda )=\frac{\mathrm{exp}\left[\frac{2A_0}{\sigma _n}\sqrt{\lambda }\mathrm{tanh}\left(\frac{\zeta \sigma _n}{2}\sqrt{\lambda }\right)\right]}{\mathrm{cosh}^{1/2}\left(\frac{\zeta \sigma _n}{2}\sqrt{\lambda }\right)}$$
(14)
and

$$G_2(\lambda )=\left(\frac{f_\lambda (\zeta )}{f_\lambda (0)}\right)^{1/2}\mathrm{exp}(\lambda ^2\sigma _n^2A_0\beta ^2\lambda A_0\zeta ^3/3),$$
(15)
where

$$\beta ^2=_0^\zeta \left[\frac{1}{f_\lambda (\zeta _1)}_{\zeta _1}^\zeta (\zeta \zeta _2)^2f_\lambda (\zeta _2)𝑑\zeta _2\right]^2𝑑\zeta _1,$$
(16)
and

$$f_\lambda (\zeta _1)=\sqrt{\zeta \zeta _1}I_{\frac{1}{4}}\left[\frac{\sqrt{\lambda }\sigma _n}{2}(\zeta \zeta _1)^2\right],$$
(17)
where $`I_\nu ()`$ is the $`\nu ^{\mathrm{th}}`$-order modified Bessel function of the first kind. The function of $`G_1(\lambda )`$ (14) is similar to the characteristic function of nonlinear phase noise Ho (2003a, b). Using $`_0^1x^{\nu +1}I_\nu (ax)𝑑x=I_{\nu +1}(a)/a`$ Gradshteyn and Ryzhik (1980), we get

$$_{\zeta _1}^\zeta (\zeta \zeta _2)^2f_\lambda (\zeta _2)𝑑\zeta _2=\frac{(\zeta \zeta _1)^{3/2}}{\sqrt{\lambda }\sigma _n}I_{\frac{3}{4}}\left(\frac{\sqrt{\lambda }\sigma _n}{2}(\zeta \zeta _1)^2\right),$$
(18)
and

$$\beta ^2=\sqrt{2}(\sqrt{\lambda }\sigma _n)^{7/2}_0^{\sqrt{\lambda }\sigma _n\zeta ^2/2}x^{1/2}I_{\frac{3}{4}}^2(x)I_{\frac{1}{4}}^2(x)𝑑x.$$
(19)
If the amplitude jitter is approximated as $`A(\zeta )A_0`$ in both variances of (4) and (5), we get

$$G_1(\lambda )\mathrm{exp}(\lambda A_0\zeta ),G_2(\lambda )\mathrm{exp}\left(\frac{1}{3}\lambda A_0\zeta ^3\right),$$
(20)
that are valid for $`A_0\sigma _n^2\zeta `$. With the approximation of (20), the characteristic functions of (10) and (11) are zero-mean Gaussian characteristic functions with variance of

$$\sigma _\mathrm{\Omega }^2(\zeta )=\frac{1}{3}A_0\sigma _n^2\zeta \text{ and }\sigma _T^2(\zeta )=\frac{1}{9}A_0\sigma _n^2\zeta ^3,$$
(21)
respectively. Note that the timing jitter has a variance increase with $`\zeta ^3`$ Gordon and Haus (1986).

From the frequency jitter of (8), the non-Gaussian distribution is induced by the term of $`\frac{1}{2}_0^\zeta w_A𝑑w_\mathrm{\Omega }`$, i.e., the noise and noise interaction. The second-order soliton perturbation also includes noise and noise interaction Kaup (1991); Haus et al. (1997). However, the equations of (1) to (3) with noise variances of (4) to (5) are directly from the first-order perturbation of soliton Kivshar and Malomed (1989); Kaup (1990); Georges (1995a). Similarly, the non-Gaussian timing jitter of (9) is induced by the term of $`\frac{1}{2}_0^\zeta (\zeta \zeta _1)w_A𝑑w_\mathrm{\Omega }`$ that also includes noise-noise interaction.

The probability density functions of frequency and Gordon-Haus timing jitter are the inverse Fourier transforms of the corresponding characteristic functions of (10) and (11), respectively. Figs. 1 plot the probability density of frequency and Gordon-Haus timing jitter as compared with the Gaussian distribution. The probability density functions are shown for signal-to-noise ratio (SNR) of $`A_0/(\sigma _n^2\zeta )=10`$ and $`20`$. The horizontal axis is normalized with respect to the standard deviation of frequency $`\sigma _\mathrm{\Omega }(\zeta )`$ and timing $`\sigma _T(\zeta )`$ jitter \[see (21)\] for Figs. 1(a) and (b), respectively. Because the characteristic functions of (10) and (11) are even real functions, the probability density functions are also even functions. Figs. 1 just plot for the positive frequency and timing jitter. Comparing the frequency jitter of Fig. 1(a) and the timing jitter of Fig. 1(b), the frequency jitter has a distribution more close to the Gaussian distribution than that of timing jitter.

Figs. 2 plot the tail probabilities corresponding to the probability density functions of Figs. 1. The tail probability is defined as $`_x^{\mathrm{}}p(x)𝑑x`$ for a probability density function of $`p(x)`$. The tail probability is compared to the complementary error function of $`\frac{1}{2}\text{erfc}\left(x/\sqrt{2}\right)`$ that is the tail probability of the Gaussian distribution in Figs. 1. Comparing the tail probability of frequency jitter of Fig. 2(a) and the timing jitter of Fig. 2(b), the frequency jitter has a distribution more close to the Gaussian distribution than that of timing jitter.

From both Figs. 1 and 2, both frequency and timing jitters have a larger tail probability than that of the Gaussian distribution. The non-Gaussian distribution leads to higher error probability than that of Gaussian distribution.

From both Figs. 1(b) and 2(b), the distribution of Gordon-Haus timing jitter has a large different to Gaussian distribution at the tail. Comparing the generating function of (15) and the Gaussian approximation of (20), the non-Gaussian timing jitter distribution is given by the factor of $`\left[f_\lambda (\zeta )/f_\lambda (0)\right]^{1/2}\mathrm{exp}(\lambda ^2\sigma _n^2A_0\beta ^2)`$ that is obviously non-Gaussian. Comparing the generating function of (14) and the Gaussian approximation of (20), the non-Gaussian frequency jitter distribution is given by two factors: the third and higher-order powers of $`\mathrm{tanh}(x)=xx^3/3+\mathrm{}`$ and the factor of $`\text{sech}^{1/2}(\zeta \sigma _n\sqrt{\lambda }/2)`$.

In conclusion, based on the first-order perturbation theory of soliton, both frequency and timing jitters are found to be non-Gaussian distributed. Amplitude, frequency, and timing jitters are all Gaussian distributed if the equations from perturbation are linearized. Without linearization, the noise projected into frequency and timing jitters are modulated by the amplitude jitter, leading to non-Gaussian distribution. The timing jitter has larger discrepancy to Gaussian distribution than that of frequency jitter.


format: text

text:
# Variational principle for frozen-in vortex structures interacting with sound waves

## I Introduction

In many physical systems the motion of continuous liquid media can be approximately described by hydrodynamic-type equations having a remarkable mathematical structure based on an underlying variational least action principle in the Lagrangian description LL2 ; ZK97 ; R99 ; R2000PRD ; R2001PRE ; RP2001PRD ; R2002 . The characteristic feature of the hydrodynamic-type systems is that they possess, besides the usual integrals of motion (the total energy, the total linear momentum, and the total angular momentum), also an infinite number of specific integrals of motion related to the freezing-in property of canonical vorticity ZK97 ; R99 ; R2001PRE . Thus, hydrodynamic equations describe an interaction between ”soft” degrees of freedom of – frozen-in vortices, and ”hard” degrees of freedom – acoustic modes. However, in the Eulerian description of flows by the density and by the velocity (or the canonical momentum) fields, vortices and sound waves are ”mixed”. Another point is that due to unresolved freezing-in constraints, the Eulerian equations of motion do not follow directly from a variational principle (see ZK97  for discussion). This work has two main purposes. The first purpose is to introduce a such general description of ideal flows, that soft and hard degrees of freedom are explicitly separated, and the frozen-in property of the vorticity is taken into account. The second purpose is to formulate a principle of least action in this representation. As the result, the acoustic waves will be described by the Eulerian fields of the fluid density and the potential component of the canonical momentum, while the canonical vorticity will be represented as a continuous distribution of frozen-in vortex lines (the so called formalism of vortex lines R99 ; R2000PRD ; R2001PRE , that previously was applied only to static density profiles). The Lagrangian of this dynamical system is a non-trivial unification of the canonical Lagrangian corresponding to purely potential flows, with a generalized Lagrangian of vortex lines.

## II Canonical formalism for fluids

### II.1 Generalized Euler equation

Typically in a complex classical system on the microscopic level there are permanently existing particles of several kinds, for instance, molecules in a gas, or the electrons and ions in a plasma. In general situation, different components can have different macroscopically averaged velocities near a same point and/or different relative concentrations in separated points. In such cases, each population of the complex fluid should be included into consideration individually, for example, as in the two-fluid plasma model discussed some later in this work. Now for simplicity we are going to consider the case when the macroscopic velocities of all components coincide and mutual relations between the concentrations are homogeneous in space and time, so the macroscopically averaged physical state of the medium at a given point $`𝒓=(x,y,z)`$ at a given time moment $`t`$ is completely determined by two quantities, namely by a scalar $`n(𝒓,t)`$, which is proportional to concentration of conservative particles of a definite sort, and by a vector $`𝒋(𝒓,t)`$, the corresponding density of flow. The field $`𝒋`$ is related by the continuity equation to the field $`n`$,

$$n_t+\text{div }𝒋=0,$$
(1)
where subscript is used to denote the partial derivative. It is clear that $`𝒋=n𝒗`$, where $`𝒗(𝒓,t)`$ is the macroscopic velocity field. Let each point of the fluid medium be marked by a label $`𝒂=(a_1,a_2,a_3)`$, so the mapping $`𝒓=𝒙(𝒂,t)`$ is the full Lagrangian description of the flow. The less exhaustive description of the flow by the fields $`n(𝒓,t)`$ and $`𝒋(𝒓,t)`$ is commonly referred as the Eulerian description. The relations between the Eulerian fields and the Lagrangian mapping are the following,

$`n(𝒓,t)={\displaystyle \delta (𝒓𝒙(𝒂,t))𝑑𝒂},`$ (2)
$`𝒋(𝒓,t)={\displaystyle \delta (𝒓𝒙(𝒂,t))𝒙_t(𝒂,t)𝑑𝒂},`$ (3)
and they satisfy the continuity equation (1) automatically.

With neglecting all dissipative processes (due to viscosity, diffusion, etc.), and assuming internal properties of the fluid homogeneous (such as the specific entropy in adiabatic flows, or the temperature in isothermal flows), the trajectories $`𝒓=𝒙(𝒂,t)`$ of fluid elements are determined by the variational principle $`\delta (\stackrel{~}{}𝑑t)/\delta 𝒙(𝒂,t)=0`$, with the Lagrangian of a special general form, actually depending only on the Eulerian fields $`n(𝒓,t)`$ and $`𝒋(𝒓,t)`$,

$$\stackrel{~}{}\{𝒙(𝒂),𝒙_t(𝒂)\}=\{n(𝒓),𝒋(𝒓)\}|_{n\{𝐱\},j\{𝐱,𝐱_t\}},$$
(4)
where the braces $`\{\mathrm{}\}`$ are used to denote functional arguments as against usual scalar or vector arguments that are denoted by the parenthesis $`(\mathrm{})`$. The equation of motion, corresponding to the Lagrangian (4), has a remarkable general structure. The usual variational Euler-Lagrange equation,

$$\frac{d}{dt}\frac{\delta \stackrel{~}{}}{\delta 𝒙_t(𝒂)}=\frac{\delta \stackrel{~}{}}{\delta 𝒙(𝒂)},$$

in the Eulerian representation has the form (generalized Euler equation)

$$\frac{}{t}\left(\frac{\delta }{\delta 𝒋}\right)=\left[\frac{𝒋}{n}\times \text{curl}\left(\frac{\delta }{\delta 𝒋}\right)\right]+\mathbf{}\left(\frac{\delta }{\delta n}\right),$$
(5)
where the variational derivative $`\delta /\delta 𝒋`$ is taken at fixed $`n(𝒓,t)`$, while the variational derivative $`\delta /\delta n`$ is taken at fixed $`𝒋(𝒓,t)`$ (compare with Refs. R99 ; R2001PRE ; R2000PRD , where this equation is written in terms of $`n`$ and $`𝒗`$ and thus looks differently). Eq.(5) together with the continuity equation (1) completely determine the time evolution of the fields $`n(𝒓,t)`$ and $`𝒋(𝒓,t)`$ .

### II.2 Hamiltonian structure

In the Hamiltonian description adopted for fluids as it is discussed in R99 ; R2001PRE , instead of the field $`𝒋`$ the variational derivative of the Lagrangian,

$$𝒑=\frac{\delta }{\delta 𝒋},$$
(6)
is used (the canonical momentum). The Hamiltonian functional is defined as the Legendre transformation,

$$\{n,𝒑\}\left(\frac{\delta }{\delta 𝒋}𝒋\right)𝑑𝒓,$$
(7)
where $`𝒋`$ should be expressed in terms of $`𝒑`$ and $`n`$. The equations of motion (5) and (1) now have the non-canonical Hamiltonian structure ZK97 ; R99 ; R2001PRE ,

$`𝒑_t=\left[{\displaystyle \frac{1}{n}}\left({\displaystyle \frac{\delta }{\delta 𝒑}}\right)\times \text{curl }𝒑\right]\mathbf{}\left({\displaystyle \frac{\delta }{\delta n}}\right),`$ (8)
$`n_t=\text{div}\left({\displaystyle \frac{\delta }{\delta 𝒑}}\right).`$ (9)
These equations can be written as $`𝒑_t=\{𝒑,\}`$ and $`n_t=\{n,\}`$, where the non-canonical Poisson bracket is given by the following expression (see ZK97 ; R99  and references therein about details),

$`\{,\}={\displaystyle \left[\frac{\delta }{\delta n}\left(\mathbf{}\frac{\delta }{\delta 𝒑}\right)\frac{\delta }{\delta n}\left(\mathbf{}\frac{\delta }{\delta 𝒑}\right)\right]𝑑𝒓}`$
$`+{\displaystyle \left(\frac{\text{curl }𝒑}{n}\left[\frac{\delta }{\delta 𝒑}\times \frac{\delta }{\delta 𝒑}\right]\right)𝑑𝒓}.`$ (10)
## III Particular examples

To emphasize universality of the employed approach, now we are going to consider several physically interesting examples.

### III.1 Eulerian hydrodynamics

Let us start with the usual Eulerian hydrodynamics. In this simple case $`n`$ is the density of the fluid, and the Lagrangian is the difference between the total macroscopic kinetic energy and the total potential energy including the thermal internal energy,

$$_E=\left(\frac{𝒋^2}{2n}\epsilon (n)nU(𝒓,t)\right)𝑑𝒓,$$

where $`\epsilon (n)`$ is the density of the internal energy, and $`U(𝒓,t)`$ is the potential of an external force.

The canonical momentum coincides with the velocity field,

$$𝒑=\frac{𝒋}{n}=𝒗,$$

and the Hamiltonian is the total energy expressed in terms of $`n`$ and $`𝒑`$,

$$_E=\left(n\frac{𝒑^2}{2}+\epsilon (n)+nU(𝒓,t)\right)𝑑𝒓.$$

The equations of motion (8-9) with this Hamiltonian take the known form

$`𝒑_t=\left[𝒑\times \text{curl }𝒑\right]\mathbf{}\left({\displaystyle \frac{𝒑^2}{2}}+\epsilon ^{}(n)+U(𝒓,t)\right),`$
$`n_t=\text{div}(n𝒑).`$
### III.2 Relativistic fluid dynamics

In the Special Relativity (with $`c=1`$ for simplicity) the field $`n`$ is the time component of the (contra-variant) 4-vector of current, while $`𝒋`$ is the space component LL6 . The absolute value $`\stackrel{~}{n}`$ of this 4-vector is equal to $`(n^2𝒋^2)^{1/2}`$ and it has the meaning of the concentration of conservative particles in the locally co-moving frame of reference. The invariant expression for the action functional implies the Lagrangian in the form (compare with R99 )

$$_r=\epsilon \left((n^2𝒋^2)^{1/2}\right)𝑑𝒓,$$

where $`\epsilon (\stackrel{~}{n})`$ is the relativistic density of the internal fluid energy including the rest energy. The canonical momentum field is defined by the relation

$$𝒑=\epsilon ^{}\left((n^2𝒋^2)^{1/2}\right)\frac{𝒋}{(n^2𝒋^2)^{1/2}}.$$

Unlike the Eulerian hydrodynamics, now it is not possible in general to get analytically the inverse relation $`𝒋(n,𝒑)`$ in order to substitute it into the expression for the Hamiltonian density,

$$h_r=\epsilon ^{}\left((n^2𝒋^2)^{1/2}\right)\frac{𝒋^2}{(n^2𝒋^2)^{1/2}}+\epsilon \left((n^2𝒋^2)^{1/2}\right).$$

Exceptions are some special dependences $`\epsilon (\stackrel{~}{n})`$ (see, for example, RP2001PRD  where the particular case $`\epsilon (\stackrel{~}{n})\stackrel{~}{n}^{4/3}`$ is considered, corresponding to the ultra-relativistic equation of state).

In the General Relativity the continuity equation is (see LL2 ; LL6 )

$$\frac{1}{\sqrt{g}}\frac{}{x^i}\left(\sqrt{g}\stackrel{~}{n}\frac{dx^i}{ds}\right)=0,$$

where $`dx^i/ds`$ is the 4-velocity of the fluid element passing through the point $`(t,𝒓)`$, and $`g=\text{det}g_{ik}`$, the determinant of the metric tensor $`g_{ik}(t,𝒓)`$. Therefore

$$n=\sqrt{g}\stackrel{~}{n}\frac{dt}{ds},j^\alpha =\sqrt{g}\stackrel{~}{n}\frac{dx^\alpha }{ds},$$

$$\stackrel{~}{n}=\sqrt{(g_{00}n^2+2g_{0\alpha }nj^\alpha +g_{\alpha \beta }j^\alpha j^\beta )}/\sqrt{g},$$
(11)
and the Lagrangian of the general-relativistic hydrodynamics is

$$_{g.r.}=\epsilon \left(\frac{\sqrt{g_{00}n^2+2g_{0\alpha }nj^\alpha +g_{\alpha \beta }j^\alpha j^\beta }}{\sqrt{g}}\right)\sqrt{g}𝑑𝒓.$$
(12)
The canonical momentum

$$p_\alpha =\epsilon ^{}(\stackrel{~}{n})\frac{(g_{0\alpha }n+g_{\alpha \beta }j^\beta )}{\sqrt{g_{00}n^2+2g_{0\alpha }nj^\alpha +g_{\alpha \beta }j^\alpha j^\beta }}$$

depends in a complicated manner on $`n`$, $`𝒋`$, $`g_{ik}`$. This circumstance makes impossible in general case to present an analytical expression for the corresponding Hamiltonian functional, but, of course, it cannot cancel the existence of the Hamiltonian in mathematical sense.

### III.3 Two-fluid plasma model

Analogously multi-component hydrodynamical models can be investigated where several fields $`n^a`$ and $`𝒋^a`$ are present corresponding to different sorts of particles, with $`a=1,2,\mathrm{},A`$. The Hamiltonian non-canonical equations of motion for such models have the same general structure as discussed above, and they should be written for each component. Below we consider a physically important example — the two-fluid plasma model. As special limit cases, this model contains the usual magnetohydrodynamics (MHD), the Electron MHD, and the Hall MHD.

#### III.3.1 Lagrangian formalism

As the start point in investigation the two-fluid plasma model, let us consider the microscopic Lagrangian of a system of electrically charged classical point particles, as it is given in the famous book by Landau and Lifshitz LL2 . This Lagrangian is approximately valid up to the second order on $`v/c`$ since excitation of the free electro-magnetic field by moving charges is negligible,

$`_{\text{micro}}={\displaystyle \underset{a}{}}{\displaystyle \frac{m_a𝒗_a^2}{2}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{ab}{}}{\displaystyle \frac{e_ae_b}{|𝒓_a𝒓_b|}}+{\displaystyle \underset{a}{}}{\displaystyle \frac{m_a𝒗_a^4}{8c^2}}`$
$`+{\displaystyle \frac{1}{4c^2}}{\displaystyle \underset{ab}{}}{\displaystyle \frac{e_ae_b}{|𝒓_a𝒓_b|}}\left(𝒗_a𝒗_b+(𝒗_a𝒏_{ab})(𝒗_b𝒏_{ab})\right),`$ (13)
where $`𝒓_a(t)`$ are the positions of the point charges $`e_a`$, $`𝒗_a(t)\dot{𝒓}_a(t)`$ are their velocities, $`𝒏_{ab}(t)`$ are the unit vectors in the direction between $`e_a`$ and $`e_b`$,

$$𝒏_{ab}=\frac{𝒓_a𝒓_b}{|𝒓_a𝒓_b|}.$$

The first double sum in Eq.(13) corresponds to the electrostatic interaction, while the second double sum describes the magnetic interaction via quasi-stationary magnetic field. It is very important that for a system with macroscopically huge number of particles the magnetic energy can be of the same order (or even larger) as the macroscopic kinetic energy produced by the first ordinary sum in Eq.(13), while the terms of the fourth order on the velocities are often negligible. Generally speaking, a large part of plasma physics is governed by this Lagrangian, at least in the cases when the velocities of particles are non-relativistic and the free electro-magnetic field is not excited significantly. Obviously, different physical problems need different procedures of macroscopic consideration of this system. The most accurate (and the most complicated) would be a kinetic description. However, for our purposes it is sufficient to apply more simple and naive procedure of the hydrodynamical averaging, that gives less accurate description of the system in terms of the concentration $`n(𝒓,t)`$ of electrons and the density $`𝒋(𝒓,t)`$ of their flow, that satisfy the continuity equation, $`n_t+\text{div }𝒋=0`$ (and corresponding fields $`N(𝒓,t)`$, $`𝑱(𝒓,t)`$ for the ions, normalized to one elementary electric charge $`e`$, so $`N=ZN_i`$.)

Neglecting all dissipative processes that take place due to collisions of the particles (though on this step we strongly reduce applicability of the following conservative two-fluid plasma model), we derive from Eq.(13) the following Lagrangian functional

$`_{2f}={\displaystyle \left[\frac{M}{2N}𝑱^2+\frac{m}{2n}𝒋^2+\frac{2\pi e^2}{c^2}[\text{curl}^1(𝑱𝒋)_{}]^2\right]𝑑𝒓}`$
$`{\displaystyle \frac{e^2}{2}}{\displaystyle \frac{d𝒓_1d𝒓_2}{|𝒓_1𝒓_2|}\left(N(𝒓_1)n(𝒓_1)\right)\left(N(𝒓_2)n(𝒓_2)\right)}`$
$`{\displaystyle \left[T_en\mathrm{ln}\frac{n}{f(T_e)}+T_i\frac{N}{Z}\mathrm{ln}\frac{N}{ZF(T_i)}\right]𝑑𝒓},`$ (14)
where $`(𝑱𝒋)_{}`$ is the divergence-free component of the total current. Here the constant $`M`$ is the ion mass per one elementary electric charge, $`M=M_i/Z`$. The electron mass is $`mM`$, and it will be neglected where possible. The magnetic energy $`(𝑩^2/8\pi )𝑑𝒓`$ is included into this Lagrangian, where the magnetic field is

$$𝑩=\frac{4\pi e}{c}\text{curl}^1(𝑱𝒋)_{}.$$
(15)
The terms with $`T_en\mathrm{ln}n`$ and $`(T_i/Z)N\mathrm{ln}N`$ (approximate expressions for the densities of the thermal free energy, LL5 ) have been introduced in order the macroscopic equations of motion to contain the pressure terms like $`p/n`$ (see the last term in Eq.(5)), where $`pnT_e`$ is the pressure of the hot electron gas, which is supposed to be isothermal with a temperature $`T_e`$. The functions $`f(T_e)`$ and $`F(T_i)`$ actually are not important since they do not contribute to the equations of motion and thus will be omitted in further equations.

It should be kept in mind that densities of the internal energy should be used instead of densities of the free energy if we suppose the flows to be isentropic. However, since the thermal conductivity is large at high temperatures ($`\kappa T_e^{5/2}/e^4m^{1/2}`$, LL10 ), the isothermal approximation usually works better than isentropic one.

#### III.3.2 Hamiltonian formalism. Hall MHD limit

For two-fluid models the Hamiltonian functional is defined as follows,

$$_{2f}\{n,𝒑,N,𝑷\}\left(\frac{\delta _{2f}}{\delta 𝒋}𝒋+\frac{\delta _{2f}}{\delta 𝑱}𝑱\right)𝑑𝒓_{2f},$$
(16)
where $`𝒋`$ and $`𝑱`$ should be expressed in terms of the electron and ion canonical momenta $`𝒑`$ and $`𝑷`$.

In our particular case we have

$`𝑷`$ $`=`$ $`{\displaystyle \frac{M𝑱}{N}}+{\displaystyle \frac{4\pi e^2}{c^2}}\text{curl}^2(𝑱𝒋)_{},`$ (17)
$`𝒑`$ $`=`$ $`{\displaystyle \frac{m𝒋}{n}}{\displaystyle \frac{4\pi e^2}{c^2}}\text{curl}^2(𝑱𝒋)_{},`$ (18)
and this results in the following approximate (valid in the limit of small $`m`$) expression for the Hamiltonian,

$`_{2f}{\displaystyle \left[\frac{N}{2M}(𝑷+𝒑)^2+\frac{n}{2m}𝒑_{}^2+\frac{c^2}{8\pi e^2}(\text{curl }𝒑)^2\right]𝑑𝒓}`$
$`+{\displaystyle \frac{e^2}{2}}{\displaystyle \frac{d𝒓_1d𝒓_2}{|𝒓_1𝒓_2|}\left(N(𝒓_1)n(𝒓_1)\right)\left(N(𝒓_2)n(𝒓_2)\right)}`$
$`+{\displaystyle [T_en\mathrm{ln}n+(T_i/Z)N\mathrm{ln}N]𝑑𝒓},`$ (19)
where $`𝒑_{}`$ is the potential component of the electron canonical momentum field.

The equations of motion for the electron component now have the non-canonical Hamiltonian structure (8-9) (and analogous equations for the ion component). Due to the large coefficient $`1/m`$ in front of $`𝒑_{}^2`$ in the Hamiltonian, $`𝒑_{}`$ should be small. A small value of $`m`$ is also the reason for the electrical quasi-neutrality on sufficiently large scales. Therefore we can put $`nN`$. Correspondingly, these simplifications will result in $`_{2f}_{simpl}\{N,𝑷+𝒑,\text{curl }𝒑\}`$. Thus, we have to deal with $`N`$, with the ion velocity

$$𝒖=(𝑷+𝒑)/M,$$

and with

$$𝝎=\text{curl }𝒑/M=(e/Mc)𝑩.$$

Let us introduce the dimensionless density

$$\rho (𝒓,t)=N(𝒓,t)/N_0,$$

the ion inertial length $`d_i`$ at the static homogeneous state,

$$d_i^2=Mc^2/(4\pi e^2N_0)=M_ic^2/(4\pi Z^2e^2N_{i0}),$$

and the speed $`s`$ of the ion-acoustic waves,

$$s^2=(T_i+ZT_e)/M_i.$$

The Hamiltonian $`_H\{\rho ,𝒖,𝝎\}=(1/M)_{simpl}`$ for this simplified two-fluid plasma model (known as the Hall MHD) is

$$_H=\left[\frac{\rho 𝒖^2}{2}+d_i^2\frac{𝝎^2}{2}+s^2\rho \mathrm{ln}\rho \right]𝑑𝒓.$$
(20)
The equations of motion are

$`_t𝒖=\left[{\displaystyle \frac{1}{\rho }}\left({\displaystyle \frac{\delta _H}{\delta 𝒖}}\right)\times \text{curl }𝒖\right]`$
$`+{\displaystyle \frac{1}{\rho }}\left[\text{curl }\left({\displaystyle \frac{\delta _H}{\delta 𝝎}}\right)\times 𝝎\right]\mathbf{}\left({\displaystyle \frac{\delta _H}{\delta \rho }}\right),`$ (21)
$`_t\rho =\text{div}\left({\displaystyle \frac{\delta _H}{\delta 𝒖}}\right),`$ (22)
$`_t𝝎=\text{curl }\left[{\displaystyle \frac{1}{\rho }}\left({\displaystyle \frac{\delta _H}{\delta 𝒖}}+\text{curl }{\displaystyle \frac{\delta _H}{\delta 𝝎}}\right)\times 𝝎\right].`$ (23)
In explicit form we easily obtain the Hall MHD equations as they are commonly known,

$$𝒖_t=\left[𝒖\times \text{curl }𝒖\right]+\frac{d_i^2}{\rho }\left[\text{curl }𝝎\times 𝝎\right]\mathbf{}\frac{𝒖^2}{2}s^2\frac{\mathbf{}\rho }{\rho },$$
(24)
$$\rho _t=\text{div}(\rho 𝒖),$$
(25)
$$𝝎_t=\text{curl}\left[\left(𝒖+\frac{d_i^2}{\rho }\text{curl }𝝎\right)\times 𝝎\right].$$
(26)
Obviously, one can normalize the length and the velocity scales to $`d_i`$ and $`s`$ respectively, that will result effectively in $`d_i=1`$, $`s=1`$.

It is interesting to note that the ion generalized vorticity can be identically equal to zero (that means $`\text{curl }𝑷=0`$, and therefore $`𝝎=\text{curl }𝒖`$), and then we have

$`𝒖_t=\left[\left(𝒖+{\displaystyle \frac{\text{curl }\text{curl }𝒖}{\rho }}\right)\times \text{curl }𝒖\right]\mathbf{}{\displaystyle \frac{𝒖^2}{2}}{\displaystyle \frac{\mathbf{}\rho }{\rho }},`$ (27)
$`\rho _t=\text{div}(\rho 𝒖).`$ (28)
This reduced system formally corresponds to a single-fluid model with the Hamiltonian $`_s\{\rho ,𝒖\}`$ as follows,

$$_s=\left[\frac{\rho 𝒖^2}{2}+\frac{(\text{curl }𝒖)^2}{2}+\rho \mathrm{ln}\rho \right]𝑑𝒓,$$
(29)
because the equations (27-28) have the standard form similar to Eqs. (8-9),

$`_t𝒖=\left[{\displaystyle \frac{1}{\rho }}\left({\displaystyle \frac{\delta _s}{\delta 𝒖}}\right)\times \text{curl }𝒖\right]\mathbf{}\left({\displaystyle \frac{\delta _s}{\delta \rho }}\right),`$ (30)
$`_t\rho =\text{div}\left({\displaystyle \frac{\delta _s}{\delta 𝒖}}\right).`$ (31)
## IV Interaction between frozen-in vortex lines and acoustic modes.

The Hamiltonian non-canonical equations (8-9) do not follow directly from a variational principle. The mathematical reason for this is a degeneracy of the corresponding non-canonical Poisson bracket (10) which is discussed, for instance, in ZK97 . The degeneracy results in frozen-in property for the canonical vorticity field $`𝝎=\text{curl }𝒑`$. However, representations of the canonical momentum in terms of auxiliary variables exist that fix topological structure of vortex lines, and then a variational formulation becomes possible. A known example of such auxiliary variables are the Clebsch variables ZK97 ; SBR2003 , when $`𝒑=\mathbf{}\phi +(\lambda /n)\mathbf{}\mu `$ and $`(n,\phi )`$, $`(\lambda ,\mu )`$ are two pairs of canonically conjugate variables. But the Clebsch representation usually is not suitable for studying localized vortex structures like vortex filaments. Below we consider another representation for the canonical momentum field, when dynamical variables are the shapes of vortex lines. For a nearly static density profile, $`nn_0(𝒓)`$, such description was used in R2001PRE ; R2000PRD  to study slow flows in spatially inhomogeneous systems. Now we are going to introduce a variational formulation valid for the general case, since the function $`n(𝒓,t)`$ is also an unknown variable. It will be demonstrated that variational principle with the Lagrangian (41) determines equations of motion for shapes of frozen-in vortex lines, for the potential component of the canonical momentum field, and for the density profile $`n(𝒓,t)`$.

So, we decompose the momentum field onto the potential component and the divergence-free component,

$$𝒑(𝒓,t)=\mathbf{}\phi (𝒓,t)+\text{curl}^1𝝎(𝒓,t).$$
(32)
Accordingly, the field $`𝒋`$ is decomposed,

$$𝒋=\frac{\delta }{\delta 𝒑}=\frac{\delta }{\delta 𝒑_{}}+\text{curl}\frac{\delta }{\delta 𝝎}𝒋_{}+𝒋_{}.$$
(33)
Obviously, the continuity equation results in the relation

$$\mathbf{}\mathrm{\Delta }^1n_t=𝒋_{}.$$
(34)
For the frozen-in vorticity field we use the so called vortex line representation. In the simplest form when the lines are closed, it reads as follows (for details and discussion see R2001PRE ),

$`𝝎(𝒓,t)={\displaystyle _𝒩}d^2\nu {\displaystyle \delta (𝒓𝑹(\nu ,\xi ,t))𝑹_\xi (\nu ,\xi ,t)𝑑\xi }`$
$`={\displaystyle \frac{𝑹_\xi (\nu ,\xi ,t)}{\text{det}𝑹/(\nu ,\xi )}}|_{𝐑=𝐫},`$ (35)
where the label $`\nu =(\nu _1,\nu _2)𝒩`$ belongs to a 2D manifold $`𝒩`$ and singles out an individual vortex line, while an arbitrary longitudinal parameter $`\xi `$ determines a point on the line. The Jacobian of the mapping $`𝑹(\nu ,\xi ,t)`$ is denoted as $`\text{det}𝑹/(\nu ,\xi )=([𝑹_{\nu _1}\times 𝑹_{\nu _2}]𝑹_\xi )`$.

The divergence-free component of the canonical momentum field now is given by the expression

$$𝒑_{}=\text{curl}^1𝝎(𝒓,t)=\frac{[𝑹_\xi \times (𝒓𝑹)]d^2\nu d\xi }{4\pi |𝒓𝑹|^3}.$$
(36)
The vorticity variation $`\delta 𝝎(𝒓,t)`$ induced by a variation $`\delta 𝑹(\nu ,\xi ,t)`$ of the vortex lines takes the form R2001PRE

$$\delta 𝝎(𝒓,t)=\text{curl}_𝐫_𝒩d^2\nu \delta (𝒓𝑹(\nu ,\xi ,t))[\delta 𝑹\times 𝑹_\xi ]𝑑\xi ,$$
(37)
which follows directly from Eq.(IV).

It should be noted that in the case of arbitrary topology of the vortex lines, one has just to replace in the above expressions $`𝑹(\nu ,\xi ,t)𝑹(𝒂,t)`$, and $`𝑹_\xi d^2\nu d\xi (𝝎_0(𝒂)\mathbf{}_a)𝑹(𝒂,t)d𝒂`$, see R2001PRE .

Eq.(37) results in the important relations R2001PRE ,

$$\frac{\delta }{\delta 𝑹}=[𝑹_\xi \times 𝒋_{}(𝑹)],$$
(38)
$$𝝎_t=\text{curl}_𝐫\left[\frac{𝑹_t\times 𝑹_\xi }{\text{det}𝑹/(\nu ,\xi )}\right]|_{𝐑=𝐫}.$$
(39)
Therefore the equation of motion for the vorticity,

$$𝝎_t=\text{curl}_𝐫[𝒗\times 𝝎],𝒗=\frac{1}{n}\frac{\delta }{\delta 𝒑},$$

means

$$\left[\frac{𝑹_t\times 𝑹_\xi }{\text{det}𝑹/(\nu ,\xi )}\right]|_{𝐑=𝐫}=\left[\frac{𝒋_{}+𝒋_{}}{n}\times 𝝎\right]+_𝐫\mathrm{\Psi }(\nu ),$$
(40)
where $`\mathrm{\Psi }(\nu _1,\nu _2)`$ is some arbitrary function of two variables. A possible choice is $`\mathrm{\Psi }=0`$, but for general purposes we will consider below $`\mathrm{\Psi }0`$.

With using Eqs.(33), (34), (38), (40), one can verify that if the quantities $`𝑹(\nu ,\xi ,t)`$, $`n(𝒓,t)`$, and $`\phi (𝒓,t)`$ obey equations of motion corresponding to the following Lagrangian,

$`_{\text{v-s}}={\displaystyle n\phi _t𝑑𝒓}\{n,\mathbf{}\phi +\text{curl}^1𝝎\{𝑹\}\}`$
$`+{\displaystyle \delta (𝒓𝑹(\nu ,\xi ,t))([𝑹_\xi \times 𝑹_t]\mathbf{}_𝐫\mathrm{\Delta }_𝐫^1n)d^2\nu 𝑑\xi 𝑑𝒓}`$
$`{\displaystyle \mathrm{\Psi }(\nu _1,\nu _2)([𝑹_{\nu _1}\times 𝑹_{\nu _2}]𝑹_\xi )n(𝑹)d^2\nu 𝑑\xi },`$ (41)
then equations (8-9) are satisfied. Indeed, the variation of $`_{\text{v-s}}`$ by $`\delta 𝑹(\nu ,\xi ,t)`$ gives the equation

$`[𝑹_\xi \times 𝑹_t]n(𝑹)[𝑹_\xi \times 𝒋_{}(𝑹)]`$
$`={\displaystyle \frac{\delta }{\delta 𝑹}}n(𝑹)\text{det}{\displaystyle \frac{𝑹}{(\nu ,\xi )}}_𝐫\mathrm{\Psi }(\nu ),`$ (42)
which is easily recognized as Eq.(40),

$`[𝑹_t\times 𝑹_\xi ]=\left[{\displaystyle \frac{𝒋_{}(𝑹)+𝒋_{}(𝑹)}{n(𝑹)}}\times 𝑹_\xi \right]`$
$`+\text{det}{\displaystyle \frac{𝑹}{(\nu ,\xi )}}(_𝐫\mathrm{\Psi })|_{𝐫=𝐑}.`$ (43)
Variation by $`\delta n(𝒓,t)`$ results in the potential component of the Eq.(8),

$$\phi _t+\mathrm{\Delta }_𝐫^1\mathbf{}_𝐫\left[\frac{𝑹_t\times 𝑹_\xi }{\text{det}𝑹/(\nu ,\xi )}\right]|_{𝐑=𝐫}=\frac{\delta }{\delta n}+\mathrm{\Psi }.$$

Finally, the variation by $`\delta \phi (𝒓,t)`$ gives the continuity equation

$$n_t+\mathbf{}𝒋_{}=0.$$

Thus, the Lagrangian (41) gives a required variational formulation for the problem of motion and interaction between localized frozen-in vortex structures (described by the mapping $`𝑹(\nu ,\xi ,t)`$) and acoustic degrees of freedom (described by the fields $`n(𝒓,t)`$ and $`\phi (𝒓,t)`$). This variational principle definitely can serve as a basis for future approximate analytical and numerical studies dealing with reduced dynamical systems where only most relevant degrees of freedom will be taken into account.

The function $`\mathrm{\Psi }(\nu _1,\nu _2)`$ can be useful to investigate nearly stationary flows, since the effective Hamiltonian,

$$\stackrel{~}{}=+\mathrm{\Psi }(\nu _1,\nu _2)([𝑹_{\nu _1}\times 𝑹_{\nu _2}]𝑹_\xi )n(𝑹)d^2\nu 𝑑\xi ,$$

has an extremum on stationary flows with the velocity field $`𝒗`$ everywhere directed along vortex surfaces. However, one should remember that existence of globally defined vortex surfaces (and thus the function $`\mathrm{\Psi }`$) is an exceptional case in the variety of 3D vector fields. In the general case one should use $`(𝝎_0(𝒂)\mathbf{}_a)𝑹(𝒂,t)d𝒂`$ instead of $`𝑹_\xi (\nu ,\xi ,t)d^2\nu d\xi `$ in the Lagrangian (41) and no function $`\mathrm{\Psi }(\nu _1,\nu _2)`$, since the labels $`\nu `$ are not defined.

Generalization of the above theory for multi-fluid models is straightforward.

As an explicit example, below is given the Hamiltonian of the Eulerian hydrodynamics in terms of $`𝑹(\nu ,\xi ,t)`$, $`n(𝒓,t)`$, and $`\phi (𝒓,t)`$,

$`_E={\displaystyle \frac{n}{2}\left(\mathbf{}\phi +\frac{[𝑹_\xi \times (𝒓𝑹)]d^2\nu d\xi }{4\pi |𝒓𝑹|^3}\right)^2𝑑𝒓}`$
$`+{\displaystyle [\epsilon (n)+nU(𝒓,t)]𝑑𝒓}.`$ (44)
### IV.1 The case of the Hall MHD

In the Hall MHD there are two frozen-in generalized vorticity fields but only one density field. It is easy to check that in this case the corresponding Lagrangian takes the following form,

$`_{\text{HMHD}}={\displaystyle \rho \phi _t𝑑𝒓}`$
$`+{\displaystyle \delta (𝒓𝑹(\nu ,\xi ,t))([𝑹_\xi \times 𝑹_t]\mathbf{}_𝐫\mathrm{\Delta }_𝐫^1\rho )d^2\nu 𝑑\xi 𝑑𝒓}`$
$`+{\displaystyle \delta (𝒓𝑿(\mu ,\eta ,t))([𝑿_\eta \times 𝑿_t]\mathbf{}_𝐫\mathrm{\Delta }_𝐫^1\rho )d^2\mu 𝑑\eta 𝑑𝒓}`$
$`_H\{\rho ,\mathbf{}\phi +\text{curl}^1𝝎\{𝑹\}+\text{curl}^1𝛀\{𝑿\},𝝎\{𝑹\}\}`$
$`{\displaystyle \mathrm{\Psi }(\nu _1,\nu _2)([𝑹_{\nu _1}\times 𝑹_{\nu _2}]𝑹_\xi )\rho (𝑹)d^2\nu 𝑑\xi }`$
$`{\displaystyle \mathrm{\Phi }(\mu _1,\mu _2)([𝑿_{\mu _1}\times 𝑿_{\mu _2}]𝑹_\eta )\rho (𝑿)d^2\mu 𝑑\eta },`$ (45)
where the vector function $`𝑹(\nu ,\xi ,t)`$ describes the frozen-in lines of the electron generalized vorticity $`𝝎`$, while $`𝑿(\mu ,\eta ,t)`$ describes the frozen-in lines of the ion generalized vorticity $`𝛀`$. The Hamiltonian $`_H\{\rho ,𝒖,𝝎\}`$ is given by Eq.(20).

### Acknowledgments

These investigations were supported by INTAS (grant No. 00292), by RFBR , by the Russian State Program of Support of the Leading Scientific Schools, and by the Science Support Foundation, Russia.


text:
# Differential Cross Sections for Higgs Boson Production at Tevatron Collider Energies

## I INTRODUCTION

Elucidation of the dynamics responsible for the breaking of electroweak symmetry is one of the primary goals in particle physics during this decade. In the current run-II of the Fermilab Tevatron, the upgraded Collider Detector Facility (CDF) and D0 detectors will search for the neutral Higgs boson ($`h`$), the vehicle of symmetry breaking in the standard model (SM), as well as the Higgs bosons of the minimal supersymmetric standard model (MSSM). The task is challenging. The production cross sections and branching fractions into channels favorable for detection are relatively small, and the backgrounds from competing processes are large Carena:2000yx .

The fully inclusive gluon-gluon fusion subprocess Wilczek:1977zn ; Graudenz:1992pv  $`gghX`$ supplies the largest cross section at Tevatron energies. Its rate is a factor four or so greater than that of the next most important subprocess, associated production $`q\overline{q}hVX`$, where $`V=W,Z`$. In the region of modest Higgs boson masses, $`m_h<135`$ GeV, the combination of gluon fusion production with the dominant decay process $`hb\overline{b}`$ should be attractive, but its promise is compromised by the overwhelming production of $`b\overline{b}`$ pairs from strong interactions background processes. For $`m_h>135`$ GeV, the channel $`hWW^{}`$ becomes dominant, and gluon fusion followed by $`hWW^{}`$ provides a potential discovery mode at the Tevatron.

Calculations of the expected differential cross sections for production of the signal and backgrounds are important for evaluation of the required measurement accuracies and detector performance. Estimations of the expected transverse transverse momentum distributions can suggest selections in this variable that should improve background rejection. In this paper, we concentrate on $`gghX`$ and discuss the behavior of the Higgs boson transverse momentum distribution in the region of small and intermediate values of $`Q_T`$. We expect our predictions to be directly pertinent in the region $`m_h>135`$ GeV, but we also entertain optimism that experimenters may find clever ways to make use of the dominant channel $`gghX`$ in the region $`m_h<135`$ GeV.

In the gluon fusion process, production occurs through triangle loops of colored (s)particles that couple to the Higgs boson and to gluons. In the SM, the most relevant contribution is from a loop of top ($`t`$) quarks. The coupling of gluons is simplified in the limit of large top quark mass $`m_t`$ Ellis:1975ap . The $`m_t\mathrm{}`$ approximation is valid to an accuracy of $`5`$% for $`m_h2m_t`$ Dawson:1990zj . Within this approach, the total cross section for $`gghX`$, is known to next-to-next-to-leading order (NNLO) accuracy Harlander:2002wh ; Ravindran:2003um . The large $`m_t`$ approximation serves well also for the transverse momentum distribution when $`m_h<2m_t`$ and the Higgs boson transverse momentum $`Q_T`$ is less than $`m_t`$ Ellis:1987xu . The next-to-leading order contributions are computed in Refs. deFlorian:1999zd ; Ravindran:2002dc ; Glosser:2002gm .

When $`Q_T`$ is comparable to $`m_h`$, there is only one hard momentum scale in the perturbative expansion of the cross section as a function of the strong coupling $`\alpha _s`$, and fixed-order computations in perturbative QCD are expected to be applicable. However, in the region $`Q_Tm_h`$, where the cross section is greatest, the coefficients of the expansion in $`\alpha _s`$ depend functionally on logarithms of the ratio of the two quantities, $`m_h`$ and $`Q_T`$, $`(\alpha _s/\pi )\mathrm{ln}^2(m_h^2/Q_T^2)`$. The relevant expansion parameter in the perturbative series is close to 1, and straightforward fixed-order perturbation theory is inapplicable for $`Q_Tm_h`$. All-orders resummation is the established method for mastering the large logarithmic coefficients of the expansion in $`\alpha _s`$ and for obtaining well-behaved cross sections at intermediate and small $`Q_T`$ Collins:1981uk ; Collins:1984kg ; Ellis:1997ii . Using renormalization group techniques, Collins, Soper, and Sterman (CSS) Collins:1984kg  devised a $`b`$ space resummation formalism that resums all logarithmic terms as singular as $`(1/Q_T^2)\mathrm{ln}^n(m_h^2/Q_T^2)`$ when $`Q_T0`$. This formalism has been used widely for computations of the transverse momentum distributions of Higgs bosons Catani:vd ; Hinchliffe:1988ap ; Kauffman:1991jt ; Yuan:1991we ; deFlorian:2000pr ; Berger:2002ut ; Bozzi:2003jy  and other processes.

In this paper, we use the $`b`$ space resummation approach discussed in some detail in our Ref. Berger:2002ut  to derive predictions for the $`Q_T`$ spectrum of Higgs boson production at the Fermilab Tevatron energy $`\sqrt{S}=1.96`$ TeV. We resum the large logarithmic terms associated with soft gluon emission to all orders in $`\alpha _s`$ obtaining well defined predictions for the full range of $`Q_T`$. We employ expressions for the parton-level hard-scattering functions valid through first-order in $`\alpha _s`$, including contributions from the glue-glue, quark-glue, and quark-antiquark incident partonic subprocesses. We present differential cross sections for masses of the Higgs boson that span the range of present interest in the SM, from $`m_h=M_Z=91.187`$ GeV to $`m_h=200`$ GeV. To illustrate interesting differences, we also provide results for the $`Z`$ boson. In this paper, we also show the integrated $`Q_T`$ distributions $`\sigma (Q_T>Q_{Tmin})`$ for rapidity within the interval $`|y|<1.0`$. This distribution indicates immediately what fraction of the Higgs boson cross section may be sacrificed if a selection is made on $`Q_T`$ for background rejection.

Our predictions are presented and discussed in Sec. II. Conclusions and a discussion of background rejection are found in Sec. III.

## II Predictions

Our predictions are based on the formalism described in detail in our Ref. Berger:2002ut . For the perturbative $`A_g`$ and $`B_g`$ functions in the Sudakov factor, we make an expansion valid through second order ($`n=2`$) in the strong coupling $`\alpha _s`$. For the short-distance coefficient functions $`C_{ab}`$ used to compute the modified parton densities, we retain the expansion through $`n=1`$, i.e., to first order in $`\alpha _s`$. We remark that the functions $`A_g`$, $`B_g`$, and $`C_{ab}`$ have well behaved perturbative expansions, free of large logarithmic terms. We use a next-to-leading order form for $`\alpha _s(\mu )`$ and next-to-leading order normal parton densities $`\varphi (x,\mu )`$ Lai:1999wy . In the fixed-order perturbative expressions that enter the “$`Y`$” function, we select a fixed renormalization/factorization scale $`\mu =\kappa \sqrt{m_h^2+Q_T^2}`$, with $`\kappa =0.5`$. For the resummed term, we take as our central value $`\mu =c/b`$, with $`c=2e^{\gamma _E}1`$. We examine sensitivity to the choice of $`\mu `$ by computing cross sections with other choices of this scale. The extrapolation into the non-perturbative region of large $`b`$ is accomplished with the form devised by Qiu and Zhang Qiu:2000hf .

In Fig. 1, we show the predicted differential $`Q_T`$ distributions. We present results for three choices of mass of the Higgs boson, $`m_h=M_Z=91.187`$ GeV, $`m_h=125`$ GeV, and $`m_h=150`$ GeV (where the $`WW^{}`$ decay channel is dominant). In all cases, the solid lines represent the predictions of resummation at next-to-leading order (NLO) accuracy and, for comparison, the dashed lines show the results of resummation at leading-order (LO) accuracy NLO-LO .

Two points are evident in the comparison of $`Z`$ boson and Higgs boson production, with $`m_h=M_Z`$, Fig. 1(c) and (d). The peak in the $`Q_T`$ distribution occurs at a smaller value of $`Q_T`$ for $`Z`$ production. At $`y=0`$, the curve peaks at $`Q_T3.0`$ GeV for $`Z`$ production and at $`Q_T7.3`$ GeV for Higgs boson production. Second, the distribution is narrower for $`Z`$ production, falling to half its maximum by $`Q_T9.0`$ GeV, whereas the half-maximum for Higgs production is not reached until $`Q_T21`$ GeV. The physics behind these important differences is that the larger QCD color factors produce more gluonic showering in the glue-glue scattering subprocess that dominates inclusive Higgs boson production than in the fermionic subprocesses relevant for $`Z`$ production. After all-orders resummation, the enhanced showering suppresses the large-$`b`$ (small $`Q_T`$) region more effectively for Higgs boson production.

For Higgs boson production, there is a substantial difference in the predictions at NLO and at LO, but the differences are slight for $`Z`$ production. The role of the $`C`$ function is responsible for these effects. The first two coefficients in the expansion of $`C`$ are published Kauffman:1991jt ; Yuan:1991we ; deFlorian:2000pr . In powers of ($`\alpha _s/\pi `$), the $`C^{(1)}`$ coefficients for the $`gg`$ and $`q\overline{q}`$ cases are

$`C_{gg}^{(1)}(z)`$ $`=`$ $`\delta (1z){\displaystyle \frac{1}{2}}\left[C_A\left({\displaystyle \frac{2\pi ^2+11}{6}}\right)+{\displaystyle \frac{\pi ^2}{6}}C_A\right].`$
$`C_{qq}^{(1)}(z)`$ $`=`$ $`\delta (1z){\displaystyle \frac{1}{2}}\left[C_F\left({\displaystyle \frac{2\pi ^224}{6}}\right)+{\displaystyle \frac{\pi ^2}{6}}C_F\right]+{\displaystyle \frac{1}{2}}C_F(1z).`$ (1)
In the $`q\overline{q}`$ case relevant for $`Z`$ production, there is a strong cancellation, $`2\pi ^224`$, reducing the effect of the NLO contribution. For Higgs boson production, the analogous term is additive, $`2\pi ^2+11`$. The constant terms $`24`$ and $`11`$, respectively, are associated with virtual corrections at NLO and can be traced to the interference between the tree and one-loop diagrams. The $`2\pi ^224`$ term is effectively the finite piece of the one-loop correction to the $`q\overline{q}Z`$ vertex while the $`2\pi ^2+11`$ term results from the finite piece of the one-loop correction to the effective $`ggh`$ vertex.

Integrating our $`Q_T`$ distributions over all $`Q_T`$ and rapidity $`y`$, we obtain the total cross sections as a function of Higgs boson mass shown in Fig. 2. We may compare these cross sections to the values listed in Table 3 of Ravindran et al Ravindran:2003um  based on their fully inclusive calculation, without resummation. Their results are obtained with different parton densities and are quoted at $`\sqrt{S}=2`$ TeV, making the comparison somewhat imprecise. The integral of our NLO resummed cross section lies between the NLO and NNLO inclusive results, slightly above the NLO values, for all masses shown, as would be expected qualitatively. Our calculation includes the ingredients $`A_g^{(1)}`$ and $`B_g^{(1)}`$ of a NLO calculation with, in addition, other ingredients such as $`A_g^{(2)}`$ of the NNLO total cross section.

Comparison of the predicted $`Q_T`$ distributions for Higgs boson production at different masses of the Higgs boson shows that the peak of the distribution shifts to greater $`Q_T`$ as $`m_h`$ grows and that the distribution broadens somewhat. At $`y=0`$, the peaks are centered at about $`Q_T=`$ 7.3, 8.2, 8.7, and 9.5 GeV for $`m_h=M_Z`$, 125, 150, and 200 GeV, respectively. The change of the $`Q_T`$ distribution with $`m_h`$ can be examined quantitatively with plots of the mean value $`<Q_T>`$ and of the root-mean-square $`<Q_T^2>^{\frac{1}{2}}`$. We find that $`<Q_T>`$ grows from about 19 GeV at $`m_h=M_Z`$ to about 28 GeV at $`m_h=200`$ GeV. The curve may be approximated with a straight line over this range, with $`<Q_T>0.079m_h+12`$ GeV. The root-mean-square of the distribution grows from about 27 GeV to about 38 GeV, reflecting the broadening of the $`Q_T`$ distribution with $`m_h`$.

For comparison with Higgs boson production, we quote our predictions for $`Z`$ production: $`<Q_T>=13`$ GeV and $`<Q_T^2>^{1/2}=20`$ GeV. The difference $`<Q_T^h><Q_T^Z>6`$ GeV at $`m_h=M_Z`$ is a manifestation of more significant gluonic radiation in Higgs boson production.

In anticipation of a discussion of selections on $`Q_T`$ to improve signal to background ratios, it is instructive to examine production of $`Z^{}`$ bosons, with mass greater than $`M_Z`$. For $`Z^{}`$ masses above $`M_Z`$ we assume that the dominant $`q\overline{q}`$ production model is unchanged except for the difference in mass of the $`Z^{}`$. We find that the change in the $`Q_T`$ distribution with mass is much less significant than is seen for the Higgs boson. For example, the peak position at $`y=0`$ increases to only 3.4 GeV at $`M_Z^{}=200`$ GeV from its value of 3 GeV at $`M_Z`$.

We introduce a new distribution for Higgs boson production that represents the integral of the differential cross section for all $`Q_T`$ greater than a minimum value, and for rapidity integrated over the interval $`|y|<1.0`$. The results are shown in Fig. 3. At $`Q_{Tmin}=10`$ GeV, the integrated distribution has dropped to 40% of its total for $`Z`$ production but to only 65%, 71%, and 74% for Higgs boson production at $`m_h=M_Z,125`$ GeV, and $`150`$ GeV, respectively. The irreducible background for Higgs boson decay to pairs of $`W`$’s arise from the QCD annihilation subprocess $`q\overline{q}WW`$. This subprocess has the same initial state structure as $`Z(Z^{})`$ boson production. The harder $`Q_T`$ spectrum for the glue-glue dominated Higgs boson production shows that the signal to background ratio can be enhanced if Higgs bosons are selected with large $`Q_T`$.

Choices of parameters are made in obtaining our results, including the renormalization/factorization scale $`\mu `$ and the non-perturbative input. Our default value $`\mu =c/b`$ with $`c=2e^{\gamma _E}1`$ provides a scale that varies inversely with the impact parameter $`b`$. This selection has the virtue that logarithmic dependence on $`\mu `$ is removed from the coefficient functions $`C`$ Collins:1984kg ; Berger:2002ut . The integration in $`b`$ space is dominated by the region near the peak in the $`b`$ distribution, $`b0.08`$ for $`m_h=125`$ GeV. Using the conjugate relationship $`Q_T1/b`$, we note that the typical hard scale $`\mu `$ is therefore $`<\mu >12.5`$ GeV ($`0.1m_h`$). To examine sensitivity to the selection of $`\mu `$, we consider other choices that are independent of $`b`$. Taking a value as large as $`\mu =2m_h`$ produces changes in the magnitude of $`d\sigma /dydQ_T`$ at the peak position of no more than about 20%, much less than the difference between the NLO and LO results in Fig. 1. Uncertainties associated with non-perturbative physics in the region of large $`b`$ are at most 1 to 2% depending on the size of the power corrections Qiu:2000hf .

## III CONCLUSIONS

Discovery of the Higgs particle is essential to shed light on the mechanism of electroweak symmetry breaking. The partonic subprocess $`g+ghX`$ dominates Higgs boson production in hadronic reactions when the Higgs boson mass is in the expected range $`m_h<200`$ GeV. The two-scale nature of the production dynamics, with mass $`m_h`$ and transverse momentum $`Q_T`$ both potentially large, and the fact that the fixed-order perturbative QCD contributions are singular as $`Q_T0`$, necessitates all-orders resummation of large logarithmic contributions in order to obtain meaningful predictions for the $`Q_T`$ distribution particularly in the region of modest $`Q_T`$ where the cross section is largest. We perform this resummation of multiple soft-gluon emission using an impact parameter $`b`$-space formalism. The resummed $`Q_T`$ distributions are determined primarily by the perturbatively calculated $`b`$-space distributions at small $`b`$, with negligible contributions from the non-perturbative region of large $`b`$.

In this paper, we present predictions for the $`Q_T`$ distributions of Higgs boson and $`Z`$ production at $`\sqrt{S}=1.96`$ TeV. At the same mass, $`m_h=M_Z`$, the predicted mean value $`<Q_T>`$ is about 6 GeV larger for Higgs boson production than for $`Z`$ boson production. For the Higgs boson, $`<Q_T>`$ grows from about 19 GeV at $`m_h=M_Z`$ to about 28 GeV at $`m_h=200`$ GeV, and the root-mean-square $`<Q_T^2>^{1/2}`$ from about 27 GeV to about 38 GeV.

Searches for the Higgs boson in its decay to $`WW`$ require an excellent understanding of the production characteristics of both the signal and backgrounds. The irreducible backgrounds arise from the QCD annihilation subprocess $`q\overline{q}WW`$. Of interest to us is the expected $`Q_T`$ dependence of the ratio of signal to irreducible background. In this paper we provide predictions for the transverse momentum distribution of the signal. The annihilation subprocess $`q\overline{q}WW`$ has the same initial state structure as $`Z(Z^{})`$ boson production. As we show, the $`Q_T`$ spectrum of $`Z^{}`$ production is predicted to be significantly softer than that for Higgs boson production. We suggest therefore that a selection of events with large $`Q_T^{WW}`$ should help in improving the signal to background ratio. The placement of the cut requires, of course, appropriate optimization to maintain signal significance. In the mass range $`m_h<135`$ GeV where $`hb\overline{b}`$ is the leading decay mode, the QCD subprocess $`ggb\overline{b}`$ supplies the dominant background. Since the signal and background are both produced from $`gg`$ scattering, we expect them to have very similar $`Q_T`$ dependences at $`m_h=m_{b\overline{b}}`$, and means other than selections on $`Q_T`$ will be necessary for background suppression in this mass interval.

## ACKNOWLEDGMENTS

Research in the High Energy Physics Division at Argonne is supported by the United States Department of Energy, Division of High Energy Physics, under Contract W-31-109-ENG-38. JWQ is supported in part by the United States Department of Energy under Grant No. DE-FG02-87ER40371.


text:
# Measurement of the Charge Asymmetry in 𝐵→𝐾^∗(892)^±𝜋^∓

Present address: \]McGill University, Montréal, Québec, Canada H3A 2T8 CLEO Collaboration

(April 22, 2003)

## Abstract

We report on a search for a $`CP`$-violating asymmetry in the charmless hadronic decay $`BK^{}(892)^\pm \pi ^{}`$, using 9.12 $`\mathrm{fb}^1`$ of integrated luminosity produced at $`\sqrt{s}=10.58`$ GeV and collected with the CLEO detector. We find $`𝒜_{CP}(BK^{}(892)^\pm \pi ^{})=0.26_{0.34}^{+0.33}`$(stat.)$`{}_{0.08}{}^{}{}_{}{}^{+0.10}`$(syst.), giving an allowed interval of $`[0.31,0.78]`$ at the 90% confidence level.

preprint: CLNS 03/1823preprint: CLEO 03-07

The Standard Model predicts that $`CP`$-violating phenomena are governed solely by the single imaginary parameter of the Cabibbo-Kobayashi-Maskawa matrix Kobayashi:fv  of complex quark couplings. The first observations of $`CP`$ violation in the neutral $`B`$ system were recently reported sin2beta , and they have been interpreted widely as induced by $`B^0`$-$`\overline{B}^0`$ mixing. To date, direct $`CP`$ violation has only been observed in the neutral kaon system directKaonCPV . Direct $`CP`$ violation in a given decay requires contributions from two or more amplitudes which differ in both $`CP`$-violating (weak) and $`CP`$-conserving (strong) phases. In the $`B`$ system, these conditions are expected to be met in some charmless hadronic decays, and direct $`CP`$ violation can occur at sizeable levels, depending on the magnitude of the strong phase difference acpTheory ; kstpiAcpTheory  or on the presence of new physics newPhysics . Previous analyses, mainly focusing on two-pseudoscalar final states, have not observed direct $`CP`$ violation in these decays cleoAcpPP ; acpExperiment . In this Report, we present a search for direct $`CP`$ violation in the vector-pseudoscalar decay $`BK^{}(892)^\pm \pi ^{}`$. We express the difference between the decay rates for $`\overline{B}^0K^{}(892)^{}\pi ^+`$ and $`B^0K^{}(892)^+\pi ^{}`$ in terms of an asymmetry, $`𝒜_{CP}`$, defined as

$$𝒜_{CP}\frac{(\overline{B}^0K^{}(892)^{}\pi ^+)(B^0K^{}(892)^+\pi ^{})}{(\overline{B}^0K^{}(892)^{}\pi ^+)+(B^0K^{}(892)^+\pi ^{})}.$$
(1)
We consider both $`K^{}(892)^\pm `$ submodes, $`K^{}(892)^\pm K_S^0\pi ^\pm `$ and $`K^{}(892)^\pm K^\pm \pi ^0`$ by analyzing the final states $`K_S^0\pi ^\pm h^{}`$ and $`K^\pm h^{}\pi ^0`$, where $`h^{}`$ denotes a charged pion or kaon. We perform a maximum likelihood fit in the $`K_S^0\pi ^\pm h^{}`$ and $`K^\pm h^{}\pi ^0`$ Dalitz plots to distinguish $`BK^{}(892)^\pm \pi ^{}`$ from other intermediate resonances or non-resonant three-body decays. The $`CP`$-averaged branching fraction for $`BK^{}(892)^\pm \pi ^{}`$ has been measured by the Belle belle  and CLEO cleo  Collaborations, and the work described in this Report is an extension of that previous CLEO analysis.

The data sample used in this analysis was produced in symmetric $`e^+e^{}`$ collisions at the Cornell Electron Storage Ring (CESR) and collected with the CLEO detector in two configurations, known as CLEO II cleonim  and CLEO II.V iivnim . It comprises 9.1 $`\mathrm{fb}^1`$ of integrated luminosity collected on the $`\mathrm{{\rm Y}}(4S)`$ resonance, corresponding to $`9.7\times 10^6`$ $`B\overline{B}`$ pairs, of which $`6.3\times 10^6`$ were taken with CLEO II.V. An additional 4.4 $`\mathrm{fb}^1`$ collected below the $`B\overline{B}`$ production threshold is used to study non-$`B\overline{B}`$ backgrounds. Of this latter luminosity, 2.8 $`\mathrm{fb}^1`$ were collected with CLEO II.V. The response of the experimental apparatus is studied with a detailed GEANT-based geant  simulation of the CLEO detector, where the simulated events are processed in a fashion similar to data.

In CLEO II, the momenta of charged particles are measured with a tracking system consisting of a six-layer straw tube chamber, a ten-layer precision drift chamber, and a 51-layer main drift chamber, all operating inside a 1.5 T superconducting solenoid. The main drift chamber also provides a measurement of specific ionization energy loss ($`dE/dx`$), which is used for particle identification. For CLEO II.V, the six-layer straw tube chamber was replaced by a three-layer double-sided silicon vertex detector, and the gas in the main drift chamber was changed from an argon-ethane to a helium-propane mixture. Photons are detected with a 7800-crystal CsI electromagnetic calorimeter, which is also inside the solenoid. Proportional chambers placed at various depths within the steel return yoke of the magnet identify muons.

Charged tracks are required to be well-measured and to satisfy criteria based on the track fit quality. They must also be consistent with coming from the interaction point in three dimensions. Pions and kaons are identified by consistency with the expected $`dE/dx`$, and tracks that are positively identified as electrons or muons are not allowed to form the $`B`$ candidate. We form $`\pi ^0`$ candidates from pairs of photons with invariant mass within 20 MeV/$`c^2`$ (approximately 2.5 standard deviations ($`\sigma `$)) of the known $`\pi ^0`$ mass. These candidates are then kinematically fitted with their masses constrained to the known $`\pi ^0`$ mass. We also require the $`\pi ^0`$ momentum to be greater than 1 GeV/$`c`$ to reduce combinatoric background from low-momentum $`\pi ^0`$ candidates. $`K_S^0`$ candidates are selected from pairs of tracks with invariant mass within 10 MeV/$`c^2`$ or $`2.5\sigma `$ of the known $`K_S^0`$ mass. In addition, $`K_S^0`$ candidates are required to originate from the beam spot and to have well-measured displaced decay vertices.

We identify $`B`$ meson candidates by their invariant mass and the total energy of their decay products. We calculate a beam-constrained mass by substituting the beam energy ($`E_b`$) for the measured $`B`$ candidate energy: $`M\sqrt{E_b^2𝐩_B^2}`$, where $`𝐩_B`$ is the $`B`$ candidate momentum. Performing this substitution improves the resolution of $`M`$ by one order of magnitude, to about 3 MeV/$`c^2`$. We define $`\mathrm{\Delta }EE_1+E_2+E_3E_b`$, where $`E_1`$, $`E_2`$, and $`E_3`$ are the energies of the $`B`$ candidate daughters. For final states with a $`K_S^0`$ and two charged tracks, the $`\mathrm{\Delta }E`$ resolution is about 20 MeV for CLEO II and 15 MeV for CLEO II.V. A $`\pi ^0`$ in the final state degrades this resolution by approximately a factor of two. $`\mathrm{\Delta }E`$ is always calculated assuming the $`h^{}`$ is a pion. Therefore, the $`\mathrm{\Delta }E`$ distribution for pions is centered at zero, while that for kaons is shifted by at least $`40`$ MeV. We accept $`B`$ candidates with $`M`$ between 5.2 and 5.3 GeV/$`c^2`$ and with $`|\mathrm{\Delta }E|`$ less than 300 MeV for $`K^\pm h^{}\pi ^0`$ and 200 MeV for $`K_S^0\pi ^\pm h^{}`$. This region includes the signal region and a high-statistics sideband for background normalization. We reject candidates that are consistent with the exclusive $`bc`$ transitions $`BD\pi `$, where $`DK\pi `$, and $`B\psi K^0`$, where $`\psi \mu ^+\mu ^{}`$ and the muons are misidentified as pions.

The main background in this analysis arises from $`e^+e^{}q\overline{q}`$, where $`q=u,d,s,c`$. To suppress this background, we calculate the angle $`\theta _{\mathrm{sph}}`$ between the sphericity axis sphericity  of the tracks and showers forming the $`B`$ candidate and that of the remainder of the event. Because of their two-jet structure, continuum $`q\overline{q}`$ events peak strongly at $`|\mathrm{cos}\theta _{\mathrm{sph}}|=1`$, while the more isotropic $`B\overline{B}`$ events are nearly flat in this variable. By requiring $`|\mathrm{cos}\theta _{\mathrm{sph}}|<0.8`$, we reject 83% of the continuum background while retaining 83% of signal $`B`$ decays. Additional separation of signal from $`q\overline{q}`$ background is provided by a Fisher discriminant Fisher  $``$ formed from eleven variables: the angle between the sphericity axis of the candidate and the beam axis, the ratio of Fox-Wolfram moments $`H_2/H_0`$ r2 , and the scalar sum of the visible momentum in nine $`10^{}`$ angular bins around the candidate sphericity axis. In the likelihood fit, we also make use of the angle between the $`B`$ candidate momentum and the beam axis, $`\theta _B`$. Angular momentum conservation causes $`B`$ mesons produced through the $`\mathrm{{\rm Y}}(4S)`$ to exhibit a $`\mathrm{sin}^2\theta _B`$ dependence, while candidates from continuum are flat in $`\mathrm{cos}\theta _B`$.

In both the $`K_S^0\pi ^\pm h^{}`$ and $`K^\pm h^{}\pi ^0`$ topologies, the $`h^{}`$ refers to the faster of the two tracks, which typically has momentum above 1 GeV/$`c`$. Because $`dE/dx`$ still provides limited separation of pions and kaons above 1 GeV/$`c`$, we make use of the $`dE/dx`$ information in the likelihood fit. In $`BK^{}(892)^\pm \pi ^{}`$ decays, this higher-momentum track is the one that recoils from the $`K^{}(892)^\pm `$, more than 99.99% of the time. Its charge uniquely distinguishes $`\overline{B}^0K^{}(892)^{}\pi ^+`$ from $`B^0K^{}(892)^+\pi ^{}`$. Thus, the charge asymmetry $`𝒜_+`$, formed using the charge of this higher-momentum track, is essentially the same as $`𝒜_{CP}`$.

Our loose selection criteria result in samples consisting primarily of background events and containing 11893 candidates for $`K_S^0\pi ^\pm h^{}`$ and 28589 for $`K^\pm h^{}\pi ^0`$. To extract yields and $`CP`$ asymmetries, we perform an unbinned maximum likelihood fit using the observables $`M`$, $`\mathrm{\Delta }E`$, $``$, $`\mathrm{cos}\theta _B`$, the two Dalitz plot variables in each topology, and the $`dE/dx`$ of the $`h^{}`$ (the faster of the two primary tracks). At high momentum, charged pions and kaons are statistically separated by their $`dE/dx`$ and by $`\mathrm{\Delta }E`$, each of which provides discrimination at the $`2.0\sigma `$ level ($`1.7\sigma `$ for CLEO II), and we fit for both $`\pi `$ and $`K`$ hypotheses simultaneously. Charged pions and kaons with momentum below 1 GeV/$`c`$ are cleanly identified by $`dE/dx`$ consistency at the $`3\sigma `$ level. The free parameters in the fit are yields ($`N`$) summed over charge states, $`N_{h^+}+N_h^{}`$, and charge asymmetries, $`𝒜_+(N_{h^+}N_h^{})/(N_{h^+}+N_h^{})`$.

The probability for a candidate to be consistent with a given component is the product of the probability density function (PDF) values for each of the input variables (neglecting correlations). The likelihood for each candidate is the sum of probabilities over the components in the fit, with relative weights determined by maximizing the total likelihood of the sample, which is given by the following expression:

$$=\underset{i=1}{\overset{\mathrm{candidates}}{}}\left\{\underset{j=1}{\overset{\mathrm{components}}{}}\left[f_j(1\pm 𝒜_+^j)\underset{k=1}{\overset{\mathrm{variables}}{}}𝒫_{ijk}\right]\right\},$$
(2)
where the $`\pm `$ refers to the charge of $`h^\pm `$ in each candidate. The $`𝒫_{ijk}`$ are the per candidate PDF values, and the $`f_j`$ and $`𝒜_+^j`$ are the free parameters optimized by the fit. The products $`f_j(1\pm 𝒜_+^j)/2`$ are constrained to sum to the fraction of candidates in the fit with the appropriate charge of $`h^\pm `$. Since the PDFs are normalized to unit integral over the fit domain, the $`f_j`$ can be interpreted as component fractions. The parameters of the $`dE/dx`$ PDFs are measured from $`DK^\pm \pi ^{}`$ decays in data. For all other variables, the signal and the background $`bc`$ PDFs are determined from high-statistics Monte Carlo samples, and the continuum PDFs are determined from data collected below the $`B\overline{B}`$ production threshold. The impact of correlations among the input variables is reduced by determining the PDFs as a function of the event location in the Dalitz plot, for coarse bins in the $`M^2(K\pi )`$-$`M^2(\pi \pi )`$ plane. We use Monte Carlo simulation to estimate the systematic uncertainty associated with neglecting any remaining correlations.

Events from $`BK_S^0\pi ^\pm h^{}`$ and $`BK^\pm h^{}\pi ^0`$, including $`BK^{}(892)^\pm \pi ^{}`$, are modeled in the fit as follows. We consider various $`B`$ decay channels with intermediate resonances ($`K^{}(892)`$, $`K_0^{}(1430)`$, $`\rho (770)`$, and $`f_0(980)`$) as well as non-resonant phase space decay. The Dalitz plot PDFs include our knowledge of the helicity structure in these decays. We neglect interference among these signal processes and assign a systematic uncertainty estimated from Monte Carlo simulation. The decays $`B𝒦^\pm h^{}`$, where $`𝒦^\pm `$ denotes $`K^{}(892)^\pm `$ or $`K_0^{}(1430)^\pm `$, are accessible through different $`𝒦^\pm `$ submodes in both the $`K_S^0\pi ^\pm h^{}`$ and $`K^\pm h^{}\pi ^0`$ topologies. To maximize our sensitivity to these decays, we fit these two topologies simultaneously, with the branching fraction and charge asymmetry for each $`𝒦^\pm h^{}`$ decay constrained to be equal in its two $`𝒦^\pm `$ submodes, which are related by isospin.

The charge symmetry of the CLEO detector, the track reconstruction software, and the $`dE/dx`$ measurement has previously been verified cleoAcpPP ; Coan:1999kh . The charge asymmetries of the fit samples are $`0.010\pm 0.009`$ for $`BK_S^0\pi ^\pm h^{}`$ and $`0.006\pm 0.006`$ for $`BK^\pm h^{}\pi ^0`$. Detection efficiencies and crossfeed among the charge-summed signal yields are measured from Monte Carlo simulated events and are accounted for in the fit for $`𝒜_+`$. We find the charge asymmetry of the detection efficiencies to be consistent with the expected null result. Crossfeed among different charge states is not included in the fit, and its effect is estimated with Monte Carlo simulation.

We perform the fit with differing combinations of intermediate resonant and non-resonant states, with up to twelve signal components. The fitted value of $`𝒜_{CP}`$ does not depend heavily on the number of signal components in the fit, and we include a systematic uncertainty for these variations. We also allow for four background components: pion and kaon hypotheses for $`h^\pm `$ for continuum background and background from $`bc`$ decays. We do not fit for charge asymmetries in the background components, but we measure them to be consistent with zero. The $`BK^{}(892)^\pm \pi ^{}`$ event yields were measured to be cleo  $`12.6_{3.9}^{+4.6}`$ for $`K^{}(892)^\pm K_S^0\pi ^\pm `$ and $`6.1_{1.9}^{+2.2}`$ for $`K^{}(892)^\pm K^\pm \pi ^0`$ with a combined statistical significance of $`4.6\sigma `$. In the fit, these yields are corrected for efficiency and crossfeed from other modes, and the $`CP`$ asymmetry in $`BK^{}(892)^\pm \pi ^{}`$ is measured to be $`𝒜_{CP}=0.26_{0.34}^{+0.33}`$$`{}_{0.08}{}^{}{}_{}{}^{+0.10}`$, where the uncertainties are statistical and systematic, respectively. The dominant contributions to the latter are statistical uncertainties in the PDFs and variations in the fitting method.

We determine the dependence of the likelihood function on $`𝒜_{CP}`$ by repeating the fit at several fixed values of $`𝒜_{CP}`$. By convoluting this function with the systematic uncertainties and integrating the resultant curve in the physical region, we construct a 90% confidence level interval of $`0.31<𝒜_{CP}<0.78`$, where the excluded regions on both sides each contain 5% of the integrated area. Figure 1 shows the likelihood function given by the fit and the effect of including systematic uncertainties.

In summary, we have measured the $`CP`$ asymmetry in $`BK^{}(892)^\pm \pi ^{}`$ using a simultaneous maximum likelihood fit to the $`BK_S^0\pi ^\pm h^{}`$ and $`BK^\pm h^{}\pi ^0`$ topologies. We obtain the value $`𝒜_{CP}=0.26_{0.34}^{+0.33}`$$`{}_{0.08}{}^{}{}_{}{}^{+0.10}`$, which is consistent with the theoretical predictions kstpiAcpTheory  of $`0.19`$ to 0.47. We also establish a 90% confidence level interval of $`[0.31,0.78]`$.

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. M. Selen thanks the Research Corporation, and A.H. Mahmood thanks the Texas Advanced Research Program. This work was supported by the National Science Foundation, and the U.S. Department of Energy.