text:
# Unified description of floppy and rigid rotating Wigner molecules formed in quantum dots

(18 November 2003)

## Abstract

Restoration of broken circular symmetry is used to explore the characteristics of the ground states and the excitation spectra of rotating Wigner molecules (RWM’s) formed in two-dimensional parabolic $`N`$-electron quantum dots. In high magnetic fields, the RWM’s are floppy rotors with the energies of the magic angular momentum $`(L)`$ states obeying $`aL+b/L^{1/2}`$. Under such fields the ground-state energies (referenced to the kinetic energy in the lowest Landau level) approach the electrostatic energy of $`N`$ point charges in the classical equilibrium molecular configuration. At zero field and strong interelectron repulsion, the RWM’s behave like quasiclassical rigid rotors whose energies vary as $`L^2`$.

The belief that the physics of two-dimensional (2D) semiconductor quantum dots (QD’s) in high magnetic fields $`(B)`$ is described jai ; haw  by composite-fermions, introduced jai2  originally for the bulk fractional quantum Hall effect (FQHE), has been recently challenged yl1  by a proposal that the fundamental physical entity in strongly correlated QD’s is a collectively rotating Wigner (or electron) molecule (RWM or REM). Indeed, it has been lately found that electrons in 2D QD’s can undergo in the strongly correlated regime a spontaneous phase transition akin to the Wigner crystallization in the bulk, forming yl2 ; yl3 ; yl4 ; mak  specific polygonal geometric structures that are called Wigner molecules (WM’s). These geometric structures break the rotational symmetry (for symmetry restoration, see below) and reflect the presence of many-body crystalline correlations that are revealed, as illustrated previously, in the “humps” of the conditional probabilities yl4 ; mak  and of the broken-symmetry electron densities yl2 ; yl3 .

The majority of recent studies address the properties of static WM’s (as a function of the QD parameters) reu ; kai ; sun ; sza . Here, we study the properties of the rotating WM. We show through microscopic many-body investigations that the RWM in high $`B`$ is a floppy rotor, while at large $`R_W`$ (and zero magnetic field) it transforms into a rigid rotor. \[The Wigner parameter $`R_WQ/\mathrm{}\omega _0`$, where $`Q`$ is the Coulomb interaction strength; $`Q=e^2/(\kappa l_0)`$, with $`l_0=\sqrt{\mathrm{}/(m^{}\omega _0)}`$ being the spatial extent of the lowest single-electron wave function in the external parabolic confinement of frequency $`\omega _0`$, and $`\kappa `$ is the dielectric constant.\] The ability to capture the physics of the electrons in QD’s in both the high magnetic field and the field-free regimes is an essential demonstration of the powerful unification offered by the RWM picture, a property not shared by other suggested approaches.

The collective rotation of the WM is inherent and natural to the molecular picture. In particular, we show that it is manifested in characteristic energy-vs-angular-momentum relations for the excitation spectra and in specific limiting values for the ground-state energies. These relations are important in the development of theories of electrons in QD’s under the influence of a magnetic-field and in the field-free case, and they can be employed as diagnostic tools for assessing the validity and applicability of alternative theoretical descriptions.

Description of the broad variation of the collective properties of electrons in 2D QD’s requires a highly flexible and accurate many-body method. In principle, exact diagonalization (EXD) could have been used; however, its computational limitations are a major obstacle. We demonstrate that the recently proposed yl3 ; yl1  two-step method of circular symmetry breaking at the unrestricted Hartree-Fock (UHF) level and of subsequent symmetry restoration via post-Hartree-Fock projection techniques is an accurate and computationally efficient approximation which provides a unified microscopic description of the emergent picture of RWM’s in 2D QD’s.

The restoration of broken circular symmetry is necessary for a proper description of RWM’s. Indeed, while the broken symmetry UHF solutions describe static WM’s reu ; kai ; sun ; sza , it is the symmetry restoration step which describes the rotation of the WM’s and underlies the differentiation between rigid and floppy (at high $`B`$) rotors.

The two-step method. In general, the localized broken symmetry orbitals are determined numerically via a selfconsistent solution of the UHF equations yl3 . An efficient alternative, however, is to approximate these orbitals by appropriate analytical expressions yl1 ; kai ; sun ; sza . Since we focus here on the second step (restoration of the circular symmetry) and the case of high $`B`$, it will be sufficient to approximate the UHF orbitals (first step of our procedure) by (parameter free) displaced Gaussian functions note5 ; namely, for an electron localized at $`Z_j`$, we use the orbital

$$u(z,Z_j)=\frac{1}{\sqrt{\pi }\lambda }\mathrm{exp}\left(\frac{|zZ_j|^2}{2\lambda ^2}i\phi (z,Z_j;B)\right),$$
(1)
with $`z=x+iy`$, $`Z_j=X_j+iY_j`$, and $`\lambda =\sqrt{\mathrm{}/m^{}\mathrm{\Omega }}`$; $`\mathrm{\Omega }=\sqrt{\omega _0^2+\omega _c^2/4}`$, where $`\omega _c=eB/(m^{}c)`$ is the cyclotron frequency. The phase guarantees gauge invariance in the presence of a perpendicular magnetic field and is given in the symmetric gauge by $`\phi (z,Z_j;B)=(xY_jyX_j)/2l_B^2`$, with $`l_B=\sqrt{\mathrm{}c/eB}`$ being the magnetic length. We only consider the case of fully polarized electrons, which is appropriate at high $`B`$.

We take the $`Z_j`$’s to coincide with the equilibrium positions (forming nested regular polygons note9 ) of $`N`$ classical point charges inside an external parabolic confinement of frequency $`\omega _0`$, and proceed to construct the UHF determinant $`\mathrm{\Psi }^{\text{UHF}}[z]`$ out of the orbitals $`u(z_i,Z_j)`$’s, $`i,j=1,\mathrm{},N`$. Correlated many-body states with good total amgular momenta $`L`$ can be extracted yl3  from the UHF determinant using projection operators; the projected energies are given by yl3

$$E_{\text{PRJ}}(L)=_0^{2\pi }h(\gamma )e^{i\gamma L}𝑑\gamma /_0^{2\pi }n(\gamma )e^{i\gamma L}𝑑\gamma ,$$
(2)
with $`h(\gamma )=\mathrm{\Psi }^{\text{UHF}}(0)|H|\mathrm{\Psi }^{\text{UHF}}(\gamma )`$ and $`n(\gamma )=\mathrm{\Psi }^{\text{UHF}}(0)|\mathrm{\Psi }^{\text{UHF}}(\gamma ),`$ where $`\mathrm{\Psi }^{\text{UHF}}(\gamma )`$ is the original UHF determinant rotated by an azimuthal angle $`\gamma `$ and $`H`$ is the many body Hamiltonian (including the vector potential, external confinement, and Coulomb two-body repulsion). We note that the UHF energies are simply given by $`E_{\text{UHF}}=h(0)/n(0)`$.

Compared to an EXD calculation, the CPU time for calculating the projected energies \[see Eq. (2)\] increases much slower as a function of $`L`$ and $`N`$. In addition, the computational efficiency of Eq. (2) is greatly enhanced by the knowledge kai  of the analytical forms of the (in general complex) matrix elements of the single-particle and two-body components of $`H`$ between the displaced Gaussians.

To test the accuracy of our method, we display in Fig. 1 the total energies $`E2\mathrm{}\mathrm{\Omega }`$ as a function of the magnetic field for the triplet state of $`N=2`$ electrons and for all three levels of calculations, i.e., the UHF, the subsequent projection \[PRJ, see Eq. (2)\], and the exact result yl4 . We observe that the UHF energies fare poorly compared to the PRJ ones; note that for $`B>7`$ T, the PRJ and exact results are practically the same. We further note that the UHF energies do not reproduce the characteristic oscillations of the exact ones; these oscillations originate from the fact that only states with magic angular momenta can become ground states (in the case of the triplet state for $`N=2`$, the magic angular momenta are $`L=2m+1`$, $`m=0,1,2,\mathrm{}`$).

Ground states at high magnetic field. Fig. 2 displays the ground state energies $`EN\mathrm{}\mathrm{\Omega }`$ for $`N=3`$ and $`N=5`$ electrons as a function of $`1/B`$. Since our method allows us to reach magnetic field values as high as $`B=400`$ T, we conclude unequivocally that, as $`B\mathrm{}`$, both the PRJ and UHF energies approach the limiting value corresponding to the classical energy of $`N`$ point-like electrons at equilibrium \[in the (0,N) configuration for $`N=25`$\] inside an external parabolic confinement of frequency $`\omega _0`$, i.e.,

$$E_{\text{cl}}^{\text{st}}(N)=(3/8)(2R_W)^{2/3}NS_N^{2/3}\mathrm{}\omega _0,$$
(3)
with $`S_N=_{j=2}^N\left(\mathrm{sin}[(j1)\pi /N]\right)^1`$ note8 . Except for the $`B\mathrm{}`$ limit, the UHF energies are higher in value than the projected ones.

The projected ground states have magic angular momenta that vary (from one trough to the other) by steps of $`N`$ units ($`N=3,4,5`$) in accordance with the trend found for a much smaller range of angular momenta in earlier EXD studies mak ; yang . With our projection method, we were able to reach remarkably high values of $`L=369,690`$, and 1100 for $`N=3,4`$, and 5, respectively. The fact that $`E_{\text{PRJ}}N\mathrm{}\mathrm{\Omega }`$ and $`E_{\text{exact}}N\mathrm{}\mathrm{\Omega }`$ (see Ref. yngl ) tend to $`E_{\text{cl}}^{\text{st}}(N)`$ for $`B\mathrm{}`$ is a significant finding that has not appeared in the FQHE literature jai ; yang ; instead it has been implied that $`EN\mathrm{}\mathrm{\Omega }0`$ as $`B\mathrm{}`$. In any case, it is doubtful that the composite-fermion approach can conform to this rigorous limit yngl ; indeed the precise value of $`E_{\text{cl}}^{\text{st}}(N)`$ depends crucially on the long-range character of the Coulomb force acting between localized classical particles, and thus it is a result innate to the WM picture.

Furthermore, this limit points out to interesting experimental ramifications. Namely, at high $`B`$, the addition energies (expressed by the finite second difference $`\mathrm{\Delta }_2E(N)`$), while not conforming to the capacitance$``$Coulomb blockade model yang , will tend to a constant value $`\mathrm{\Delta }_{2,\text{cl}}^{\text{st}}E(N)=E_{\text{cl}}^{\text{st}}(N+1)+E_{\text{cl}}^{\text{st}}(N1)2E_{\text{cl}}^{\text{st}}(N)`$ (e.g., 3.446 meV for $`N=4`$ and the parameters listed in the caption of Fig. 2); this constant value maybe detected experimentally.

The yrast rotational band at high B. The yrast band consists of the lowest energy states (for a given $`B`$) having magic angular momenta (the ground state is a member of this band and the remaining states are excited ones). The yrast band for $`N=5`$ calculated at $`B=60`$ T is displayed in Fig. 3.

The energy of the yrast states is composed jai  of two contributions; namely, $`E(L)=E^c(L)+V(L)`$, where $`E^c(L)`$ is the confinement energy associated with the external potential (including the kinetic energy), and $`V(L)`$ is the interaction energy due to the Coulomb force. For high $`B`$, it is well known jai ; yang  that the confinement energy varies linearly with $`L`$; in particular, $`E^c(L)\mathrm{}(\mathrm{\Omega }\omega _c/2)L+N\mathrm{}\mathrm{\Omega }`$.

We have calculated the projected values \[see Eq. (2)\] for these two components, namely, $`E_{\text{PRJ}}^c(L)`$ and $`V_{\text{PRJ}}(L)`$, for $`N=26`$ and a variety of $`B`$ values and QD parameters. Our ability to calculate efficiently high $`L`$ values allowed us to confirm numerically the linear $`L`$-dependence of the confinement energy $`E^c(L)`$, and also the precise proportionality coefficient $`\mathrm{}(\mathrm{\Omega }\omega _c/2)`$. Moreover, rather unexpectedly, we also find that the interaction energy behaves asymptotically as $`L^{1/2}`$, i.e., $`V(L)C_V/L^{1/2}`$, where $`C_V`$ is a constant given below.

TABLE I lists $`V_{\text{PRJ}}(L)`$ for $`N=5`$ and $`B=60`$ T, as well as for $`N=6`$ \[the (1,5) configuration\] and $`B\mathrm{}`$; the latter limit can be easily taken by setting yl1  $`\lambda =l_B\sqrt{2}`$ in Eq (1), which restricts the electrons to the lowest Landau level. From TABLE I, one can see that the quantity $`fV_{\text{PRJ}}(L)L^{1/2}/C_V`$ approaches unity as the angular momentum intreases, where $`C_V=K^{3/2}(S_K/4+\delta _{K,N1})e^2/(\kappa \lambda )`$, with $`K=N`$ or $`N1`$ for WM’s in the $`(0,N)`$ or $`(1,N1)`$ configurations, respectively.

The $`C_V/L^{1/2}`$ asymptotic dependence of the interaction energy can be readily associated with a model mak  of a classical floppy molecule rotating inside a parabolic confinement characterized by a frequency $`\mathrm{\Omega }`$ note2 . Indeed, the energy of such a molecule is given by

$`E_{\text{cl}}^{\text{rot}}(K)=`$ $`\mathrm{}^2L^2/(2𝒥(a))+𝒥(a)\mathrm{\Omega }^2/2\mathrm{}\omega _cL/2`$ (4)
$`+K(S_K/4+\delta _{K,N1})e^2/(\kappa a),`$
where $`𝒥(a)=Km^{}a^2`$ is the moment of inertia of $`K`$ point-like electrons located at the vertices of a regular polygon of radius $`a`$. The second term is the potential energy due to the confinement, and the last term is the classical Coulomb-repulsion electrostatic energy. At given $`B`$, the radius of this floppy molecule varies with $`L`$, and for large angular momentum (and/or high $`B`$) it is given by $`a\lambda \sqrt{L/K}`$ \[$`a\sqrt{2\mathrm{}cL/(eBK)}`$ for $`B\mathrm{}`$\]. Substitution into Eq. (4) yields the aforementioned linear and $`1/L^{1/2}`$ contributions, i.e.,

$$E_{\text{cl}}^{\text{rot}}(K)\mathrm{}(\mathrm{\Omega }\omega _c/2)L+C_V/L^{1/2}.$$
(5)
A semiclassical approximation $`L_{\text{gs}}^{\text{cl}}`$ of the ground-state angular momentum minimizes Eq. (5), yielding $`L_{\text{gs}}^{\text{cl}}B`$. Consequently, the ground-state radius $`a_{\text{gs}}r_0`$ as $`B\mathrm{}`$, where $`r_0=l_0R_W^{1/3}(S_K/4+\delta _{K,N1})^{1/3}`$ is the radius of the static classical Wigner molecule. Note that, in the ground state as $`B`$ increases, the WM rotates faster in order to maintain the requisite value of $`a_{\text{gs}}=r_0`$.

We note that this inverse-square-root-of-$`L`$ law at constant $`B`$ for the Coulomb interaction energy has also been overlooked in the FQHE literature jai ; yang . This law provides a nontrivial test in that it lends further support to the RWM picture; especially since it also applies in the $`B\mathrm{}`$ limit (see right part of TABLE I), which restricts the electrons to the lowest Landau level (FQHE regime). We note that high $`L`$’s correspond to lower fractional fillings, experimentally achieveable in high-mobility samples.

The rigid rotor at zero magnetic field. At $`B=0`$, we found it advantageous to allow the width $`\lambda `$ and the positions $`Z_j`$’s of the displaced Gaussians to vary in order to minimize (for each $`L`$) the projected energy \[Eq. (2)\]. For large $`R_W`$’s and for $`N=25`$ electrons, we find that the energies of the states in the yrast band can be approximated by

$$E_{\text{PRJ}}(L)E_{\text{PRJ}}(0)+C_RL^2,$$
(6)
where the rigid-rotor coefficient $`C_R`$ is essentially a constant whose value is very close (see the $`\stackrel{~}{f}`$ values) to the classical one, corresponding to point charges in their $`(0,N)`$ equilibrium configuration inside a parabolic confinement of frequency $`\omega _0`$; i.e., $`C_RC_R^{\text{cl}}=\mathrm{}^2/(2𝒥(r_0))`$.

In Table 2, we list calculated $`E_{\text{PRJ}}(L)`$ values for $`N=5`$ when $`R_W=200`$. The fact that $`\stackrel{~}{f}C_R/C_R^{\text{cl}}1`$ for $`L45`$ illustrates that the RWM behaves as a quasiclassical rigid rotor. For smaller values of $`R_W`$, the rigidity of the RWM is progressively reduced yl4 . We note that $`E_{\text{PRJ}}(0)N\mathrm{}\omega _0`$ is very close to the classical electrostatic value given by Eq. (3).

Conclusions. Using the method of restoration of broken circular symmetry, we have shown that the RWM’s in QD’s exhibit characteristic properties as a result of their collective rotation: (I) in high $`B`$, they behave like floppy rotors and the energies of the magic states in the lowest rotational (yrast) band, associated with magic angular momenta, obey $`aL+b/L^{1/2}`$. In addition, the ground-state energies (referenced to the kinetic energy in the lowest Landau level) approach the limit of the electrostatic energy of $`N`$ point charges in their classical equilibrium molecular configuration; (II) at $`B=0`$ and strong interelectron repulsion, the RWM’s behave like quasiclassical rigid rotors whose yrast band exhibits an $`L^2`$-dependence. The rotating-Wigner-molecule picture unifies the description of the various physical regimes of strongly correlated electrons in 2D QD’s. This unification is achieved through the computationally powerful method of restoration of broken symmetries via projection techniques.

This research is supported by the U.S. D.O.E. (Grant No. FG05-86ER45234).


format: text

text:
# Multidimensional solitons in a low-dimensional periodic potential

## Abstract

Using the variational approximation(VA) and direct simulations, we find stable 2D and 3D solitons in the self-attractive Gross-Pitaevskii equation (GPE) with a potential which is uniform in one direction ($`z`$) and periodic in the others (but the quasi-1D potentials cannot stabilize 3D solitons). The family of solitons includes single- and multi-peaked ones. The results apply to Bose-Einstein condensates (BECs) in optical lattices (OLs), and to spatial or spatiotemporal solitons in layered optical media. This is the first prediction of mobile 2D and 3D solitons in BECs, as they keep mobility along $`z`$. Head-on collisions of in-phase solitons lead to their fusion into a collapsing pulse. Solitons colliding in adjacent OL-induced channels may form a bound state (BS), which then relaxes to a stable asymmetric form. An initially unstable soliton splits into a three-soliton BS. Localized states in the self-repulsive GPE with the low-dimensional OL are found too.

PACS numbers: 03.75.Fi, 05.30.Jp, 05.45.-a

Introduction. Solitons in multidimensional nonlinear Schrödinger equations (NLSEs) with a periodic potential have recently attracted considerable interest. In particular, self-trapping of spatial beams in nonlinear photonic crystals is described by a 2D equation. In this case, despite the possibility of collapse , simulations reveal robust 2D solitons in the self-focusing model . A similar medium can be created by a grid of laser beams illuminating a photorefractive sample .

Similar 2D and 3D models with a periodic potential describe a Bose-Einstein condensate (BEC) trapped in an optical lattice (OL) ; this application is relevant, as experimental techniques for loading BECs into multidimensional OLs were recently developed . Stable solitons can be supported by an OL even in self-repulsive BECs . In the case of self-attraction, 2D and 3D solitons (including 2D vortices) are stable in the self-focusing model with the OL potential , despite the possibility of the collapse.

An issue of direct physical relevance, and a subject of the present work, are multidimensional solitons in media equipped with periodic potentials of a lower dimension, viz., quasi-1D (Q1D) and Q2D lattices in the 2D and 3D cases, respectively. In optics, the 2D equation in the spatial domain governs the beam propagation in a layered bulk medium along the layers, which is an extension of a 1D multichannel system introduced in Ref. , where the potential was induced by transverse modulation of the refractive index (RI). In the temporal domain, the 2D and 3D equations govern, respectively, the longitudinal propagation of spatiotemporal optical solitons in a layered planar waveguide, or in a bulk medium with the RI periodically modulated in both transverse directions. The 2D and 3D cases directly apply to BECs loaded in a Q1D or Q2D lattice. The physical significance of this setting is two-fold: first, in the experiment it is much easier to create a lower-dimensional lattice than a full-dimensional one, both in BECs and in optics, hence the use of such lattices is the most straightforward way to create multidimensional solitons in these media. Second, the solitons created this way can freely move in the unconfined direction, which also suggests a possibility to study their collisions, and to look for their bound states (BSs). As yet, no other way to create multidimensional mobile solitons in BECs and their BSs has been proposed. In optics, solitons of this type suggest new applications. Indeed, in an optical medium with the full-dimensional periodic potential, transfer of a trapped beam from one position to another is difficult, as the necessary external push strongly disturbs the beam . In the lower-dimensional potential, the beam can slide along the free direction, making the transfer easy. In BECs confined by a low-dimensional OL, matter-wave solitons can be driven in the free direction by a weak laser beam.

The model and VA. The normalized form of the self-focusing 2D NLSE with a Q1D periodic potential of the strength $`\epsilon `$ is

$$iu_t+u_{xx}+u_{yy}+[\epsilon \mathrm{cos}(2x)+\chi |u|^2]u=0,$$
(1)
where $`\chi =\pm 1`$ in the case of the self-attraction/repulsion. In BECs, $`t`$ is time, while in optics it is the propagation distance. For BECs or spatial optical solitons, $`x`$ and $`y`$ are transverse coordinates; for spatiotemporal optical solitons in a 2D waveguide with anomalous chromatic dispersion, $`y`$ is a temporal variable. The 3D version of Eq. (1) is

$$iu_t+u_{xx}+u_{yy}+u_{zz}+\left\{\epsilon \left[\mathrm{cos}(2x)+\mathrm{cos}(2y)\right]+\chi |u|^2\right\}u=0.$$
(2)
Equations (1) and (2) conserve the Hamiltonian, the norm $`N=\left|u(𝐫)\right|^2𝑑𝐫`$ (the number of atoms in the BEC, or total power/energy of the spatial/spatiotemporal optical soliton), and the momentum along the free direction. The equations are normalized so that the period of the potential is $`\pi `$, the free parameters being $`\epsilon `$ and $`N`$.

Stationary solutions to Eq. (1) are looked for as $`u(x,y,t)=U(x,y)\mathrm{exp}(i\mu t)`$, with a chemical potential $`\mu `$ (propagation constant, in optics), which yields

$$\mu U+U_{xx}+U_{yy}+\left[\epsilon \mathrm{cos}(2x)+\chi U^2\right]U=0,$$
(3)
and similar in the 3D case. To apply the VA (variational approximation) to Eq. (3), we adopt the ansatz $`U=A\mathrm{exp}\left[\left(ax^2+by^2\right)/2\right]`$, or $`U=A\mathrm{exp}\left[\left(a(x^2+y^2)+bz^2\right)/2\right]`$ in the 3D case, with the norms $`N_{2\mathrm{D}}=\pi A^2/\sqrt{ab}`$, and $`N_{3\mathrm{D}}=\pi ^{3/2}A^2/(a\sqrt{b})`$. Following the procedure of the application of VA to multidimensional solitons , we derive variational equations from the Lagrangian of Eq. (3). In the 2D case, they are ($`\mathrm{}e^{1/a}`$)

$$N\chi =\left(4\pi /a\right)\sqrt{a^22\epsilon \mathrm{}},\mu =a\epsilon \mathrm{}(13/a),$$
(4)
and in the 3D case \[in both cases, $`b=\left(a^22\epsilon \mathrm{}\right)/a`$\],

$$N\chi =2\left(2\pi /a\right)^{3/2}\sqrt{a^22\epsilon \mathrm{}},\mu =a/2\epsilon \mathrm{}\left(23/a\right).$$
(5)
Evidently, solutions are possible only if $`\chi >0`$ (self-attraction). The 2D solutions exist in the interval

$$N_{\mathrm{min}}4\pi \sqrt{1\epsilon /\epsilon _{\mathrm{cr}}}<N<N_{\mathrm{max}}4\pi $$
(6)
\[see Fig. 1(a)\], if $`\epsilon <\epsilon _{\mathrm{cr}}e^2/80.924`$ ($`N_{\mathrm{min}}=0`$ if $`\epsilon >\epsilon _{\mathrm{cr}}`$), while 3D solutions can be found for any $`N`$. The fixed lattice period rules out establishing general scaling relations for these soliton families; in the case of very large $`\epsilon `$, the potential valleys become isolated, then one is actually dealing with cigar-shaped  Q1D solitons, that obey the usual 1D NLSE scaling.

The solution families which definitely meet the condition $`d\mu /dN<0`$ are expected to be stable pursuant to the Vakhitov-Kolokolov (VK) criterion , while ones with $`d\mu /dN`$ positive or nearly zero should be unstable. It follows from Eqs. (4) and (5) and is evident in Figs. 1(b,c) that the 2D solitons are stable in the existence interval (6). At small $`\epsilon `$, the 3D solutions are unstable for any $`N`$. Equation (5) predicts that a stability interval of a width $`\mathrm{\Delta }N\sqrt{\epsilon \epsilon _0}`$ appears around $`N_010.1\pi ^{3/2}`$ if $`\epsilon `$ exceeds $`\epsilon _00.242`$.

Numerical results. Solitons are generated by the imaginary-time evolution method , starting with the VA-predicted waveforms. Stability is verified by direct simulations of perturbed solitons in real time. Typical examples of the thus found stable 2D and 3D solitons are displayed in Fig. 2 and 3, the 3D ones being nearly isotropic in the $`(x,y)`$ plane and elongated in the free direction $`z`$. As is seen in Fig. 1(a), the existence limits for the stable 2D solitons are close to the VA prediction, unless $`N`$ is small; in that case, the discrepancy is due to the fact that the soliton develops a multi-peaked structure \[Fig. 2(b)\], which is not accounted for by the simple ansatz adopted above. The ansatz can be generalized for this case, setting $`U=\left(A+B\mathrm{cos}x\right)\mathrm{exp}\left[\left(ax^2+by^2\right)/2\right]`$, but the final result is then quite messy. A numerically found stability region for the 3D solitons is also quite close to the VA prediction, and multi-peaked 3D solitons relate to their single-peaked counterparts (Fig. 3) the same way as Fig. 2(b) to Fig. 2(a).

The VA predicts that, in the 3D case with the Q1D (rather than Q2D) lattice, all the solitons are VK-unstable. Accordingly, simulations never produced stable solitons in this case. This feature can be explained by the fact that, in the free 2D subspace, the soliton is essentially tantamount to the unstable Townes soliton . In the 2D NLSE with the 2D lattice, stable solitons with intrinsic vorticity were also found in Ref. . Vortices were found in the present 2D model too, but they are always unstable.

We have also investigated the case of self-repulsion, $`\chi =1`$, when the low-dimensional lattice potential cannot support a completely localized pulse. However, adding a usual parabolic trap readily gives rise to stable states, which are solitons across the lattice and Thomas-Fermi states along the free direction. It is noteworthy that, with $`\chi <0`$, these states assume a multi-peaked or single-peaked shape under the action of a weaker or stronger OL potential, respectively, which is opposite to what was reported above for $`\chi >0`$, cf. Fig. 2. An explanation is that, in the case of the self-repulsion and $`\epsilon =0`$, no solitons exist (even unstable ones, like the above-mentioned Townes soliton).

The influence of the parabolic trap was checked too in the case of $`\chi =+1`$. If the 2D or 3D soliton is displaced from the central position, it performs harmonic oscillations along the free direction, completely preserving its integrity. Actually, the mobility in the free direction is the most essential difference of the present multidimensional solitons from those predicted in other BEC models . This suggests new possibilities, such as collisions between moving solitons.

Simulations demonstrate that 2D and 3D in-phase solitons which collide head-on with velocities $`\pm v`$, which are below a critical value $`v_\mathrm{c}`$, merge into a single pulse, whose norm exceeds the critical value $`N_{\mathrm{max}}`$ \[see Eq. (6)\], hence it quickly collapses. If $`v>v_\mathrm{c}`$, the solitons pass through each other (for instance, $`v_\mathrm{c}=8.5`$ for $`\epsilon =2.0`$, when the norm of each 2D soliton is $`N=2.5\pi `$), which is explained by the fact that the collision time, $`1/v`$, is then smaller than the collapse time, $`1/N`$. Colliding $`\pi `$-out-of-phase solitons always bounce back, and two such solitons, placed inside the trap, perform stable oscillations with periodic elastic collisions.

In the low-dimensional potential, collisions are also possible between solitons moving in adjacent “tracks” (channels). In that case, they pass each other quasi-elastically if the collision velocity is large. If it is not too large, each soliton captures a small “satellite” in the other channel. If the velocity falls below a critical value for given $`\epsilon `$, the outcome of the collision is altogether different: two solitons form a quiescent symmetric bound state (BS) \[Fig. 4(a)\]. Later, the symmetric BS develops an intrinsic instability, spontaneously rearranging itself into a completely stable asymmetric BS \[Fig. 4(b)\]. This behavior resembles the symmetry-breaking instability of optical solitons in dual-core fibers  and the macroscopic quantum-self-trapping effect in BECs confined to a double-well potential . In fact, this is the first example of a stable BS of two solitons in BECs. Formation of more complex stable 2D and 3D BSs also occurs in a different way: VK-unstable waveforms predicted by the VA undergo violent evolution, shedding $`50\%`$ of their norm and eventually forming a BS of three solitons, with a narrow tall one in the middle, and small-amplitude broad satellites in the adjacent channels \[see Fig. 5\].

Conclusion. We have demonstrated that the periodic potential whose dimension is smaller by $`1`$ than the full spatial dimension readily supports stable single- and multi-peaked solitons in 2D and 3D self-focusing media (although the quasi-1D potential cannot stabilize 3D solitons), which suggest new settings for experimental search of solitons. 2D solitons exist in a finite interval of the values of the norm, and 3D ones are stable if the lattice strength exceeds the minimum value. In the case of self-repulsion with a parabolic trap, stable localized states are found too. The dependence of their structure on the lattice strength is opposite to that in the case of the self-attraction: with the increase of the strength, a single-peaked state is changed by a multi-peaked one.

These solitons are the first example of  mobile multi-dimensional pulses in BECs. Head-on collisions between the in-phase solitons may lead to their fusion and collapse, while out-of-phase solitons collide elastically indefinitely many times. Collision between solitons in adjacent tracks may create a bound state (BS), which then relaxes to an asymmetric stable shape. Three-soliton BSs are created as a result of the evolution of initially unstable solitons.

Acknowledgements. We appreciate discussions with Y.S. Kivshar. B.B.B. thanks the Physics Department at the University of Salerno (Italy) for a research grant. B.A.M. appreciates hospitality of the same Department, and partial financial support through the Excellence-Research-Center grant No. 8006/03 from the Israel Science Foundation. M.S. acknowledges partial financial support from the MIUR, through the inter-university project PRIN-2000, and from the European LOCNET grant HPRN-CT-1999-00163.


text:
# On the Complex Wilson Coefficients of New Physics Scenarios

## I Introduction and Conclusion

At present for certain decays there are a number of impressive theoretical predictions in great agreement with the experimental results. Due to the consistency, uncertainties around one $`\sigma `$ level seems acceptable. In general, the difference between Standard Model (SM) results and experimental values are seen on the half of New Physics (NP) scenarios and used to put constraints on the new parameters(i.e. see Ref.Ref00 ). If SM predictions and experimental measurements are very close then naturally NP is orphaned. The aim of the present work is to draw attention on the side effects of bounding Wilson coefficients with the mentioned procedure.

One of the best examples of such decays is the $`BX_s\gamma `$ decay, which has solid theoretical and experimental background, widely used to test and constrain a number of NP scenarios. Recent review of this decay can be found in Ref.Ref01 . Taking into consideration NP effects, the present status of the decay related with branching ratio can be expressed as

$`\eta \pm \sigma =\alpha |C|^2,`$ (1)
where $`\eta ,\sigma `$ denotes experimental measurements and uncertainties, $`\alpha `$ satisfies proportionality related with the decay. With the definition, $`C=C^{SM}(\mu )+C^{NP}(\mu )`$, $`C^{SM,NP}(\mu )`$ are the Wilson coefficients at the low energy scale $`\mu `$.

From the experimental side, average of experimental measurements of the branching ratio for the mentioned decay can be written as $`Br(BX_s\gamma )=(\eta \pm \sigma )\times 10^4`$ with $`\eta =3.34`$, $`\sigma =0.38`$ Refs.Ref02 ; Ref03 ; Ref04 ; Ref05 ; Ref06 ; Ref07 . From the viewpoint of the SM the situation is very similar, roughly $`\eta 3.33.7`$,6 $`\sigma 0.3`$ Refs.Ref08 ; Ref09 ; Ref10 ; Ref11 ; Ref12 .

It can be dreamed that one day theoretical and experimental results are very reliable and uncertainties are under control, even $`\sigma _{Theory,Experiment}0`$ is approached. Meanwhile, $`\sigma `$ is commented as a nice room for a number of new physics scenarios. Since, it can be used by NP, to fill the gap between SM theory and experiment. In the case when $`\sigma `$ is converging to zero and SM predictions are right in the correct place, $`C^{NP}0`$ and hence speculative existence of new physics is ruled out.

However, what should be stressed is that, this criteria is not enough to reject the possibility of a hidden NP. To make it clear assume the worst situation for NP, $`C^{SM}`$ satisfies the equality given in Eq.1 without any theoretical or experimental error, whence no room for NP scenarios,

$$\frac{\eta \pm 0}{\alpha }=|C^{SM}(\mu )|^2=|C^{SM}(\mu )+C^{NP}(\mu )|^2.$$
(2)
Notice that if $`C^{NP}(\mu )`$ is a real quantity, then $`0`$ is the first solution for NP Wilson coefficient, second solution is $`2C^{SM}(\mu )`$. Once it is assumed complex in the form

$$C^{NP}(\mu )=C_r(\mu )+\mathrm{𝑖𝐶}_i(\mu ),$$
(3)
at first sight $`C_r(\mu )=C^{SM}(\mu ),C_i(\mu )=\pm C^{SM}(\mu )`$ is the solution, and the general solution for the NP scenario can be written as

$$C_r^2(\mu )+C_i^2(\mu )=2C_r(\mu )C^{SM}(\mu ),$$
(4)
lying on the complex $`C^{NP}(\mu )`$ plane. Complex phases in the final form of Wilson coefficients is possible in certain extensions of the SM. As an example let us mention one of the most popular extensions of the SM, Two Higgs Doublet Model 2HDM(III) Ref.Ref13 ; Ref14 ; Ref15 ; Ref16 . In this model, not only at the matching scale $`\mu _W`$ but also for the $`\mu _b`$ scale there are complex couplings which can be written in the appropriate form as follows

$$C_7^{eff}(\mu )=C_7^{SM}(\mu )+C_{7,r}^{2HDM}(\mu )+\mathrm{𝑖𝐶}_{\mathit{7},i}^{\mathit{2}HDM}(\mu ).$$
(5)
Now, to see the effect of this approach, for simplicity, assume theoretical predictions of the SM and experimental results are given as follows

$`Br(BX_s\gamma )_{SM}`$ $`=`$ $`(\eta _{SM}\pm 0)\times 10^4,`$
$`Br(BX_s\gamma )_{Exp}`$ $`=`$ $`(\eta _{Exp}\pm \sigma _{Exp})\times 10^4;`$
$`\eta _{Exp,SM}=3.5,\sigma _{Exp}=0.4.`$
In this approximation, Wilson coefficient of the SM takes $`C_7^{SM}\pm \{0.32,0.34,0.36\}`$ values in accordance with the error $`\sigma _{Exp}`$ is given in Eq.I. and the proportionality constant $`\alpha 30.4\times 10^4`$, is defined in Eq.1. Notice that we assume no uncertainty for $`\sigma _{SM}`$ and set $`\eta _{SM}=3.5`$, hence neglect theoretical errors on purpose. Considering new physics, allowed ranges of the new Wilson coefficients $`C_{7,r}^{NP},C_{7,i}^{NP}`$ can be extracted from Fig.1. In the figure the choice $`\sigma =0`$ is presented in the middle of both left and right regions, stressing the possible solutions of $`C_7^{NP}`$ for $`Br(BX_s\gamma )_{SM}=Br(BX_s\gamma )_{Experiment}`$. By looking from one dimension if we set $`C_{7,i}^{NP}=0`$, possible values of $`C_{7,r}^{NP}`$ can be extracted from the same figure 1. As it can be deduced from the figure, considering complex phases enriches phenomenology.

To summarize, scenarios permitting complex components in the final form of evolved Wilson coefficients have a rich potential for NP effects. Free parameters of NP should not be accepted as real from the beginning. Even if SM and experiment are in complete agreement, complex parts can help NP to survive, when we consider the issue as an alternative solution. While observing the possibility of such a structure of $`NewPhysics`$, branching ratio can not be used, solely, to refuse or accept a new model. Nevertheless, since we have well motivated theoretical and experimental background, it is possible to back-transform Eq.4. and obtain the most general form of the NP scenarios at the matching scale. The price we have to pay is, at least, two more unknowns which should be fixed by other measurements or probably best by CP asymmetry $`(A_{CP})`$ of the related decays.


text:
# Parity of Anti-Decuplet Baryons Revisited from Chiral Soliton Models

## Abstract

We recalculate masses and widths of anti-decuplet baryons in the case of positive parity from chiral soliton models, provided that the member $`\mathrm{\Xi }_{3/2}`$ of the anti-decuplet has a mass 1.86 GeV as reported recently. Calculations show that there are no convincing candidates for the nonexotic members of the anti-decuplet available in the baryon listings. Up to the leading order of $`m_s`$ and 1/$`N_c`$, the width formula for the decay of the anti-decuplet baryons to the octet depends only on $`SU(3)`$ symmetry model-independently, except the coupling constant. Similarly we give a width formula for the decay of negative parity baryons belong to certain $`SU(3)`$ baryon multiplet by pure symmetry consideration. By this formula, we find that if we have an anti-decuplet with negative parity and that the masses are the same as those given by chiral soltion models, the identification of $`N(1650)`$ as $`N_{\overline{10}}`$ are inconsistent with experiments for $`N(1650)N\pi `$ while the widths agree with other two decay channels involving strangeness. And $`\mathrm{\Sigma }(1750)`$ seems to be a reasonable candidate for $`\mathrm{\Sigma }_{\overline{10}}`$.

preprint: hep-ph/0311331

Due to recent reports LEPS ; DIAN ; CLAS ; SAPH ; HERMES  for the existence of pentaquark state $`\mathrm{\Theta }^+`$, Skyrme’s old idea Skyr  of identifying baryons as solitons has aroused unprecedent interest. In the large $`N_c`$ limit witt , the soliton picture of baryons can be proved consistent with QCD. The quantization of the SU(3) Skyrmion in collective coordinates Guad  not only gives the correct baryon octet and decuplet, but also predicts the next lighter baryon multiplet, the anti-decuplet Mano-Chem , in which $`\mathrm{\Theta }^+`$ is the lightest member. In the quark language, $`\mathrm{\Theta }^+`$ is of the minimal five-quark configuration $`|uudd\overline{s}`$ qqqqqbar ; GM99 , and thus has the exotic strangeness number S=+1. Predictions about the mass and width of $`\mathrm{\Theta }^+`$ from chiral soliton models Pra ; penta1 ; Diak ; Weig  have played an important role in the searches of $`\mathrm{\Theta }^+`$. However, there are controversies about the parity of $`\mathrm{\Theta }^+`$. Naively, the lowest $`qqqq\overline{q}`$ state should have negative parity in quark models qqqqqbar ; QMDS . Analysis from QCD sum rule QCDSR  as well as earlier lattice QCD calculation LatticeQCD  show that $`\mathrm{\Theta }^+`$ can be consistently identified as a pentaquark state provided its $`J^P=\frac{1}{2}^{}`$. However, analysis of the stability of pentaquarks $`uudd\overline{s}`$ shows that Stan  $`p`$-shell with positive parity is lower than $`s`$-shell with negative parity for pentaquarks. The quark cluster models DQM  also predict a positive parity for $`\mathrm{\Theta }^+`$. In literatures Diak ; Weig , the anti-decuplet baryons are predicted to have positive parity from chiral soliton models. However, there are also opinions anti-sol  that the success predictions from chiral soliton models could be fortuitous. Up to now, the parity of $`\mathrm{\Theta }^+`$ has not been experimentally determined. Combining the observed pentaquark $`\mathrm{\Theta }^+`$ mass from available experiments LEPS ; DIAN ; CLAS ; SAPH ; HERMES  and the pentaquark $`\mathrm{\Xi }_{3/2}`$ mass from a recent experiment NA49 , Diakonov and Petrov Diak2  re-identified N$`_{\overline{10}}`$ and $`\mathrm{\Sigma }_{\overline{10}}`$. However, there are no such two corresponding baryons with $`J^P=\frac{1}{2}^+`$ in available baryon listings PDG . Thus, they suggested that there must be two new baryons as the missing members of the anti-decuplet. And in AAPSW , the modified PWA analysis also indicates that the $`N(1710)`$ is not the appropriate candidate to be a member of the anti-decuplet, instead, $`N(1680)`$ or $`N(1730)`$ with positive parity was suggested.

The main purpose of this note is to discuss the probability of picturing $`\mathrm{\Theta }^+`$, $`\mathrm{\Xi }(1862)`$, $`N(1650)`$ and $`\mathrm{\Sigma }(1750)`$ as the members of an anti-decuplet with negative parity model-independently by pure symmetry consideration. In soliton picture, there also exist pentaquark states with one heavy anti-quark ($`\overline{Q}qqqq`$) and with $`J^P=\frac{1}{2}^{}`$ in the bound state approach RS ; OPM1 . However, in the collective quantization approach, how to describe both positive and negative parity baryons is still not solved.

In the $`SU(3)`$ chiral soliton model, the fundamental object is the chiral field $`U(x)`$

$$U(x)=\mathrm{exp}\left[i\frac{\lambda _b\varphi _b(x)}{f_\pi }\right],$$
(1)
where $`f_\pi `$93 MeV is the observed pion decay constant, $`\lambda _b`$ are the eight Gell-Mann $`SU(3)`$ matrices and $`\varphi _b(x)`$ are the eight pseudoscalar meson fields. Under the space inversion transformation, $`\varphi _b(x)`$ transforms as

$$\widehat{P}\varphi _b(𝐱,t)\widehat{P}^{}=\varphi _b(𝐱,t),$$
(2)
where $`\widehat{P}`$ is the parity operator. Accordingly, $`U(x)`$ transforms as

$$\widehat{P}U(x)\widehat{P}^{}=\mathrm{exp}\left[i\frac{\lambda _b\varphi _b(𝐱,t)}{f_\pi }\right]=U^{}(𝐱,t).$$
(3)
In chiral limit, the action of Skyrme model is of the form

$`I={\displaystyle \frac{f_\pi ^2}{4}}{\displaystyle d^4x\text{Tr}(_\mu U_\mu U^{})}`$
$`+{\displaystyle \frac{1}{32e^2}}{\displaystyle d^4x\text{Tr}[_\mu UU^{},_\nu UU^{}]^2}+N_c\mathrm{\Gamma },`$ (4)
which has a hedgehog solitonic solution under the assumption of maximal symmetry Skyr ; Guad

$$U_1(𝐱)=\left(\begin{array}{cc}\mathrm{exp}[i(\widehat{𝐫}\tau )F(r)]& \begin{array}{c}0\\ 0\end{array}\\ \begin{array}{cc}0& 0\end{array}& 1\end{array}\right),$$
(5)
where $`e`$ is introduced to stabilize the solitons by Skyrme; $`\mathrm{\Gamma }`$ is the Wess-Zumino term; $`F(r)`$ is the spherical-symmetric profile of the soliton, the solution of the nonlinear equation of motion; $`\tau `$ are the three Pauli matrices; and $`\widehat{𝐫}`$ is the unit vector in space. The action is invariant under $`SU(3)_L\times SU(3)_R`$ transformation. However, we are only interested in those $`U(x)`$ with the same vacuum at $`r\mathrm{}`$

$$U(x)=A(t)U_1(𝐱)A(t)^1,ASU(3).$$
(6)
Inserting (6) into (3) gives

$$\widehat{P}U(x)\widehat{P}^{}=U(x),$$
(7)
and

$`\widehat{P}A(t)\widehat{P}^{}=A(t),`$ (8)
$`\widehat{P}F(r)\widehat{P}^{}=F(r).`$ (9)
Then quantize the system about this solitonic solution for collective coordinates $`A`$. To leading order only symmetry modes (collective coordinates) are important, thus we only treat collective coordinates quantum mechanically, and the SU(3) symmetric effective action in the large $`N_c`$ limit leads to the collective Hamiltonian Guad :

$$\widehat{H}=M_{cl}+\frac{1}{2I_2}\left[\widehat{C}^{(2)}\frac{1}{12}(N_cB)^2\right]+\left(\frac{1}{2I_1}\frac{1}{2I_2}\right)\widehat{𝐉}^2,$$
(10)
where $`M_{cl}`$, $`I_1`$ and $`I_2`$ are given by the 3-dimensional space coordinate integrals of even functions of $`F(r)`$ and $`e`$, and are treated model-independently and fixed by experimental data in our work. $`M_{cl}`$ is the classical soliton mass; $`I_1`$ and $`I_2`$ are moments of inertia; $`\widehat{C}^{(2)}`$=$`\underset{a=1}{\overset{8}{}}\widehat{G}_a^2`$ is the quadratic (Casimir) operator of the vectorial group SU(3)v, and in the representation $`(p,q)`$, its eigenvalue $`C^{(2)}=\frac{1}{3}[p^2+q^2+pq+3(p+q)]`$; $`\widehat{G}_a(A)(a=18)`$ are the generators of SU(3)v; and $`\widehat{J}_i(A)(i=13)`$ are the the generators of the spin group SU(2)s. Using Eq. (8), we have

$`\widehat{P}\widehat{H}(A)\widehat{P}^{}=\widehat{H}(A),`$ (11)
$`\widehat{P}\widehat{G}_a(A)\widehat{P}^{}=\widehat{G}_a(A),\widehat{P}\widehat{J}(A)\widehat{P}^{}=\widehat{J}(A).`$ (12)
The wave function $`\mathrm{\Psi }_{\nu \nu ^{}}^{(\mu )}`$ of baryon $`B`$ in the collective coordinates is of the form

$$\mathrm{\Psi }_{\nu \nu ^{}}^{(\mu )}(A)=\sqrt{\text{dim}(\mu )}D_{\nu \nu ^{}}^{(\mu )}(A),$$
(13)
where $`(\mu )`$ denotes an irreducible representation of the SU(3) group; $`\nu `$ and $`\nu ^{}`$ denote $`(Y,I,I_3)`$ and $`(1,J,J_3)`$ quantum numbers collectively; $`Y`$ is the hypercharge of $`B`$; $`I`$ and $`I_3`$ are the isospin and its third component of $`B`$ respectively; $`J_3`$ is the third component of spin $`J`$; and $`D_{\nu \nu ^{}}^{(\mu )}(A)`$ are representation matrices. Eq. (8) means that the wave function (13) has positive parity. However, in the procedure of collective coordinate quantization, only collective parts are quantized, and modes orthogonal to the symmetry modes are treated classically. In our case, these modes should be related with the coordinate $`r`$. However, by Eq. (9), we know that a function of $`r`$ may still possess negative parity, this may be suggestive to extend soliton picture to negative parity baryons because the part we treat classically may still contribute negative parity. The symmetry breaking Hamiltonian Bolt  is

$$H^{}=\alpha D_{88}^{(8)}+\beta Y+\frac{\gamma }{\sqrt{3}}\underset{i=1}{\overset{3}{}}D_{8i}^{(8)}J^i,$$
(14)
where the coefficients $`\alpha `$, $`\beta `$, $`\gamma `$ are proportional to the strange quark mass and model dependent, and are treated model-independently and fixed by experiments; and $`D_{ma}^{(8)}(A)`$ is the adjoint representation of the $`SU(3)`$ group and defined as

$$D_{ma}^{(8)}(A)=\frac{1}{2}\text{Tr}(A^{}\lambda ^mA\lambda ^a).$$
(15)
We can use perturbation theory to calculate the baryon states in collective coordinates on the basis of flavor symmetry states (13) Park

$`|N=|N;8+C_{\overline{10}}|N;\overline{10}+C_{27}|N;27,`$ (16)
$`|\mathrm{\Sigma }=|\mathrm{\Sigma };8+C_{\overline{10}}|\mathrm{\Sigma };\overline{10}+{\displaystyle \frac{\sqrt{6}}{3}}C_{27}|\mathrm{\Sigma };27,`$ (17)
$`|\mathrm{\Xi }=|\mathrm{\Xi };8+C_{27}|\mathrm{\Xi };27,`$ (18)
$`|\mathrm{\Lambda }=|\mathrm{\Lambda };8+{\displaystyle \frac{\sqrt{6}}{2}}C_{27}|\mathrm{\Lambda };27,`$ (19)
$`|\mathrm{\Theta }^+=|\mathrm{\Theta }^+;\overline{10},`$ (20)
$`|N_{\overline{10}}=|N;\overline{10}C_{\overline{10}}|N;8+{\displaystyle \frac{\sqrt{30}}{80}}C_{27}^{(\overline{10})}|N;27,`$ (21)
$`|\mathrm{\Sigma }_{\overline{10}}=|\mathrm{\Sigma };\overline{10}C_{\overline{10}}|\mathrm{\Sigma };8+{\displaystyle \frac{1}{4\sqrt{5}}}C_{27}^{(\overline{10})}|\mathrm{\Sigma };27,`$ (22)
$`|\mathrm{\Xi }_{3/2}=|\mathrm{\Xi }_{3/2};\overline{10}+{\displaystyle \frac{\sqrt{6}}{16}}C_{27}^{(\overline{10})}|\mathrm{\Xi }_{3/2};27,`$ (23)
and the coefficients for other multiplets are

$`C_{\overline{10}}={\displaystyle \frac{1}{3\sqrt{5}}}(\alpha +{\displaystyle \frac{\gamma }{2}})I_2,`$ (24)
$`C_{27}={\displaystyle \frac{\sqrt{6}}{25}}(\alpha {\displaystyle \frac{\gamma }{6}})I_2,`$ (25)
$`C_{27}^{(\overline{10})}=(\alpha {\displaystyle \frac{7}{6}}\gamma )I_2,`$ (26)
where $`|B;\mu `$ denotes a flavor symmetry state $`\mathrm{\Psi }_{\nu \nu ^{}}^{(\mu )}(A)`$ with $`\nu \nu ^{}`$ denoting the quantum numbers of baryon $`B`$. From the analysis above, we see that this mixing is due to $`H^{}`$, and is irrelevant to parity. Thus this mixing is still valid if the anti-decuplet baryons have negative parity. We know that in a baryon multiplet, different baryon functions are related by $`\widehat{T}_\pm (=\widehat{G_1}\pm i\widehat{G}_2)`$, $`\widehat{V}_\pm (=\widehat{G_4}i\widehat{G}_5)`$ and $`\widehat{U}_\pm (=\widehat{G_6}\pm i\widehat{G}_7)`$, which commute with $`\widehat{P}`$ by (12). Therefore, to fix the parity, we can only find a candidate with definite parity from experiments, which has the same mass and width to the corresponding member of the anti-decuplet, and then fix the parity of other baryons in the multiplet accordingly.

Ref. NA49  first reported evidence for the existence of a narrow $`\mathrm{\Xi }^{}\pi ^{}`$ baryon resonance with mass of $`1.862\pm 0.003`$ GeV and width below the detector resolution of about $`0.018`$ GeV, and this state is considered as a candidate for the pentaquark $`\mathrm{\Xi }_{\frac{3}{2}}^{}`$ in the anti-decuplet predicted from chiral soliton models Diak . Provided $`\mathrm{\Theta }^+`$ with mass of 1.54 GeV and $`\mathrm{\Xi }_{\frac{3}{2}}`$ with mass of 1.86 GeV, we recalculate the coefficients in (10) and (14) as well as the masses of the other two members of the anti-decuplet. The results are as follows

$$\begin{array}{ccc}1/I_1=154MeV;\hfill & 1/I_2=399\text{MeV};\hfill & \mathrm{\Sigma }78MeV;\hfill \\ \alpha =663\text{MeV};\hfill & \beta =12\text{MeV};\hfill & \gamma =185\text{MeV};\hfill \\ m_{N_{\overline{10}}}=1.65\text{GeV};\hfill & m_{\mathrm{\Sigma }_{\overline{10}}}=1.75\text{GeV},\hfill & \end{array}$$
(27)
where $`\mathrm{\Sigma }`$ is the pion-nucleon Sigma term. From available baryon listings PDG , there is only one particle, $`N(1650)`$ with negative parity, for the resonances of $`N`$ with mass around 1.65 GeV and spin 1/2, and there exists $`\mathrm{\Sigma }(1750)`$ with $`J^P=\frac{1}{2}^{}`$ for the resonances of $`\mathrm{\Sigma }`$ with mass around 1.75 GeV. It has been also found wu\_ma1\_2  surprisingly that all the nonexotic members of the 27-plet with spin 3/2, calculated with the above parameters, can be identified from the available baryon listings PDG , agreeable with experiments both in mass and width.

Inspired by the above observation, we want to know whether there exists some symmetry between those four baryons of $`\mathrm{\Theta }^+`$, $`N(1650)`$, $`\mathrm{\Sigma }(1750)`$, and $`\mathrm{\Xi }(1862)`$, assuming them as a set of members in an anti-decuplet. Since we still do not know how to investigate it from chiral soliton picture, we want to check it by symmetry consideration model-independently. For the decay of $`J^{}J^++0^{}`$, we have a following effective interaction langrangian

$$_I=ig\overline{\psi }\psi \phi _m,$$
(28)
and from it we have

$`\mathrm{\Gamma }(BB^{}m)={\displaystyle \frac{G_s^2}{4\pi }}{\displaystyle \frac{m_B^{}}{m_B}}|𝐩|\left\{\begin{array}{c}\frac{\text{dim}(\mu ^{})}{\text{dim}(\mu )}\left|\begin{array}{c}\underset{\gamma }{}\left(\begin{array}{cc}8& \mu ^{}\\ Y_mI_m& Y_\rho I_\rho \end{array}|\begin{array}{c}\mu _\gamma \\ Y_\nu I_\nu \end{array}\right)\left(\begin{array}{cc}8& \mu ^{}\\ 00& 1J_\rho \end{array}|\begin{array}{c}\mu _\gamma \\ 1J_\nu \end{array}\right)\end{array}\right|^2\end{array}\right\},`$ (39)
(40)
where $`\phi _m`$ is any pseudoscalar meson, $`\overline{\psi }`$ is a baryons belonging to baryon multiplet $`(\mu )`$ with negtive parity, and $`\psi `$ is a baryons belonging to baryon multiplet $`(\mu ^{})`$ with positive parity. In this paper, we will only deal with an anti-decuplet with negative parity. Thus we get the widths of anti-decuplet baryons in Table 1 by using this formula and identifying $`\mathrm{\Theta }^+`$, $`N(1650)`$, $`\mathrm{\Sigma }(1750)`$, and $`\mathrm{\Xi }(1862)`$ as the members of the anti-decuplet, and the results calculated in the case of positive parity from chiral soliton models are still list in Table 1.

| Table 1. The widths of baryons in the anti-decuplet | | | |
| --- | --- | --- | --- |
| Mode | estimation(MeV) | $`J^P=\frac{1}{2}^{}`$(MeV) | $`J^P=\frac{1}{2}^+`$(MeV) |
| $`\mathrm{\Theta }^+KN`$ | $`<25`$ | 25$`(\text{input})`$ | 25$`(\text{input})`$ |
| N(1650)$`N\pi `$ | $`80`$171 | 12 | 25 |
| N(1650)$`N\eta `$ | $`419`$ | 8 | 8 |
| N(1650)$`\mathrm{\Lambda }K`$ | $`421`$ | 4 | 1.5 |
| $`\mathrm{\Sigma }(1750)N\overline{K}`$ | 6$``$64 | 7 | 3 |
| $`\mathrm{\Sigma }(1750)\mathrm{\Sigma }\pi `$ | $`<`$12.8 | 8 | 14 |
| $`\mathrm{\Sigma }(1750)\mathrm{\Sigma }\eta `$ | 9$``$88 | 2 | 0.2 |
| $`\mathrm{\Sigma }(1750)\mathrm{\Lambda }\pi `$ | seen | 12 | 28 |
| $`\mathrm{\Xi }_{3/2}\mathrm{\Xi }\pi `$ | $`<`$18(?) | 24 | 46 |
| $`\mathrm{\Xi }_{3/2}\mathrm{\Sigma }\overline{K}`$ |  | 18 | 15 |

In summary, we recalculated the masses and widths of the baryons in the anti-decuplet. The calculated results show that if we accept $`\mathrm{\Theta }^+`$ and the reported particle $`\mathrm{\Xi }(1862)`$ in NA49  as members of the anti-decuplet, we will get the mass of $`N_{\overline{10}}`$ around 1.65 GeV and the mass of $`\mathrm{\Sigma }_{\overline{10}}`$ around 1.75 GeV. From available baryon listings PDG , we find that there are only $`N(1650)`$ and $`\mathrm{\Sigma }(1750)`$ with negative parity. And inspired by this, we give a width formula for the negative parity baryon decay on symmetry consideration and give all the widths of baryons of an anti-decuplet with negative parity but the same masses as calculated from chiral soliton models (Table 1). From Table 1, The width of $`\mathrm{\Sigma }(1750)`$ in the case of $`J^P=\frac{1}{2}^{}`$ fits experimental data better, but for $`N(1650)`$, there is an obvious deviation in the process $`N(1650)N\pi `$ while the widths agree with other two decay channels $`N(1650)N\eta `$ and $`N(1650)\mathrm{\Lambda }K`$. There could be a possibility for a missing negative parity $`N_{\overline{10}}`$ with mass around 1650 MeV and a narrow width or that the SU(3) breaking plays a role. The width of $`\mathrm{\Xi }_{3/2}`$ with negative parity also agrees better with experimental observation. In the quark picture of the anti-decuplet, $`N_{\overline{10}}`$ and $`\mathrm{\Sigma }_{\overline{10}}`$ contain hidden quark-antiquark pairs, while $`N(1650)`$ and $`\mathrm{\Sigma }(1750)`$ are the orbital angular excitations in the quark model. Thus, to determine the parity of the anti-decuplet, we need further experiments to find missing non-exotic members and also to determine their parities.

We are grateful for helpful discussions with Yajun Mao, Zhi-Yong Zhao, Han-Qing Zheng, and Shi-Lin Zhu. This work is partially supported by National Natural Science Foundation of China under Grant Numbers 10025523 and 90103007.


text:
# Hierarchical Clustering Based on Mutual Information

## 1 Introduction

Classification or organizing of data is very important in all scientific disciplines. It is one of the most fundamental mechanism of understanding and learning \[Jain and Dubes, 1988\]. Depending on the problem, classification can be exclusive or overlapping, supervised or unsupervised. In the following we will be interested only in exclusive unsupervised classification. This type of classification is usually called clustering or cluster analysis.

An instance of a clustering problem consist of a set of objects and a set of properties (called characteristic vector) for each object. The goal of clustering is the separation of objects into groups using only the characteristic vectors. Indeed, in general only certain aspects of the characteristic vectors will be relevant, and extracting these relevant features is one field where mutual information (MI) plays a major role \[Tishby et al., 1999\], but we shall not deal with this here. Cluster analysis organizes data either as a single grouping of individuals into non-overlapping clusters or as a hierarchy of nested partitions. The first approach is called partitional clustering (PC), the second one is hierarchical clustering (HC). One of the main features of HC methods is the visual impact of the dendrogram which enables one to see how objects are being merged into clusters. From any HC one can obtain a PC by restricting oneself to a “horizontal” cut through the dendrogram, while one cannot go in the other direction and obtain a full hierarchy from a single PC. Because of their wide spread of applications, there are a large variety of different clustering methods in use \[Jain and Dubes, 1988\].

The crucial point of all clustering algorithms is the choice of a proximity measure. This is obtained from the characteristic vectors and can be either an indicator for similarity (i.e. large for similar and small for dissimilar objects), or dissimilarity. In the latter case it is convenient but not obligatory if it satisfies the standard axioms of a metric (positivity, symmetry, and triangle inequality). A matrix of all pairwise proximities is called proximity matrix. Among HC methods one should distinguish between those where one uses the characteristic vectors only at the first level of the hierarchy and derives the proximities between clusters from the proximities of their constituents, and methods where the proximities are calculated each time from their characteristic vectors. The latter strategy (which is used also in the present paper) allows of course for more flexibility but might also be computationally more costly.

Quite generally, the “objects” to be clustered can be either single (finite) patterns (e.g. DNA sequences) or random variables, i.e. probability distributions. In the latter case the data are usually supplied in form of a statistical sample, and one of the simplest and most widely used similarity measures is the linear (Pearson) correlation coefficient. But this is not sensitive to nonlinear dependencies which do not manifest themselves in the covariance and can thus miss important features. This is in contrast to mutual information (MI) which is also singled out by its information theoretic background \[Cover and Thomas, 1991\]. Indeed, MI is zero only if the two random variables are strictly independent.

Another important feature of MI is that it has also an “algorithmic” cousin, defined within algorithmic (Kolmogorov) information theory \[Li and Vitanyi, 1997\] which measures the similarity between individual objects. For a thorough discussion of distance measures based on algorithmic MI and for their application to clustering, see \[Li et al., 2001, Li et al., 2002\].

Another feature of MI which is essential for the present application is its grouping property: The MI between three objects (distributions) $`X,Y,`$ and $`Z`$ is equal to the sum of the MI between $`X`$ and $`Y`$, plus the MI between $`Z`$ and the combined object (joint distribution) $`(XY)`$,

$$I(X,Y,Z)=I(X,Y)+I((X,Y),Z).$$
(1)
Within Shannon information theory this is an exact theorem (see below), while it is true in the algorithmic version up to the usual logarithmic correction terms \[Li and Vitanyi, 1997\]. Since $`X,Y,`$ and $`Z`$ can be themselves composite, Eq.(1) can be used recursively for a cluster decomposition of MI. This motivates the main idea of our clustering method: instead of using e.g. centers of masses in order to treat clusters like individual objects in an approximative way, we treat them exactly like individual objects when using MI as proximity measure.

More precisely, we propose the following scheme for clustering $`n`$ objects with MIC:
(1) Compute a proximity matrix based on pairwise mutual informations; assign $`n`$ clusters such that each cluster contains exactly one object;
(2) find the two closest clusters $`i`$ and $`j`$;
(3) create a new cluster $`(ij)`$ by combining $`i`$ and $`j`$;
(4) delete the lines/columns with indices $`i`$ and $`j`$ from the proximity matrix, and add one line/column containing the proximities between cluster $`(ij)`$ and all other clusters;
(5) if the number of clusters is still $`>2`$, goto (2); else join the two clusters and stop.

In the next section we shall review the pertinent properties of MI, both in the Shannon and in the algorithmic version. This is applied in Sec. 3 to construct a phylogenetic tree using mitochondrial DNA and in Sec. 4 to cluster the output channels of an independent component analysis (ICA) of an electrocardiogram (ECG) of a pregnant woman, and to reconstruct from this the maternal and fetal ECGs. We finish with our conclusions in Sec. 6.

## 2 Mutual Information

### 2.1 Shannon Theory

Assume that one has two random variables $`X`$ and $`Y`$. If they are discrete, we write $`p_i(X)=\mathrm{prob}(X=x_i)`$, $`p_i(Y)=\mathrm{prob}(Y=x_i)`$, and $`p_{ij}=\mathrm{prob}(X=x_i,Y=y_i)`$ for the marginal and joint distribution. Otherwise (and if they have finite densities) we denote the densities by $`\mu _X(x),\mu _Y(y)`$, and $`\mu (x,y)`$. Entropies are defined for the discrete case as usual by $`H(X)=_ip_i(X)\mathrm{log}p_i(X)`$, $`H(Y)=_ip_i(Y)\mathrm{log}p_i(Y)`$, and $`H(X,Y)=_{i,j}p_{ij}\mathrm{log}p_{ij}`$. Conditional entropies are defined as $`H(X|Y)=H(X,Y)H(Y)=_{i,j}p_{ij}\mathrm{log}p_{i|j}`$. The base of the logarithm determines the units in which information is measured. In particular, taking base two leads to information measured in bits. In the following, we always will use natural logarithms. The MI between $`X`$ and $`Y`$ is finally defined as

$`I(X,Y)`$ $`=`$ $`H(X)+H(Y)H(X,Y)`$ (2)
$`=`$ $`{\displaystyle \underset{i,j}{}}p_{ij}\mathrm{log}{\displaystyle \frac{p_{ij}}{p_i(X)p_i(Y)}}.`$
It can be shown to be non-negative, and is zero only when $`X`$ and $`Y`$ are strictly independent. For $`n`$ random variables $`X_1,X_2\mathrm{}X_n`$, the MI is defined as

$$I(X_1,\mathrm{},X_n)=\underset{k=1}{\overset{n}{}}H(X_k)H(X_1,\mathrm{},X_n).$$
(3)
This quantity is often referred to as (generalized) redundancy, in order to distinguish it from different “mutual informations” which are constructed analogously to higher order cumulants, but we shall not follow this usage. Eq.(1) can be checked easily, together with its generalization to arbitrary groupings. It means that MI can be decomposed into hierarchical levels. By iterating it, one can decompose $`I(X_1\mathrm{}X_n)`$ for any $`n>2`$ and for any partitioning of the set $`(X_1\mathrm{}X_n)`$ into the MIs between elements within one cluster and MIs between clusters.

For continuous variables one first introduces some binning (‘coarse-graining’), and applies the above to the binned variables. If $`x`$ is a vector with dimension $`m`$ and each bin has Lebesgue measure $`\mathrm{\Delta }`$, then $`p_i(X)\mu _X(x)\mathrm{\Delta }^m`$ with $`x`$ chosen suitably in bin $`i`$, and 111Notice that we have here assumed that densities really exists. If not e.g. if $`X`$ lives on a fractal set), then $`m`$ is to be replaced by the Hausdorff dimension of the measure $`\mu `$.

$$H_{\mathrm{bin}}(X)\stackrel{~}{H}(X)m\mathrm{log}\mathrm{\Delta }$$
(4)
where the differential entropy is given by

$$\stackrel{~}{H}(X)=𝑑x\mu _X(x)\mathrm{log}\mu _X(x).$$
(5)
Notice that $`H_{\mathrm{bin}}(X)`$ is a true (average) information and is thus non-negative, but $`\stackrel{~}{H}(X)`$ is not an information and can be negative. Also, $`\stackrel{~}{H}(X)`$ is not invariant under homeomorphisms $`x\varphi (x)`$.

Joint entropies, conditional entropies, and MI are defined as above, with sums replaced by integrals. Like $`\stackrel{~}{H}(X)`$, joint and conditional entropies are neither positive (semi-)definite nor invariant. But MI, defined as

$$I(X,Y)=𝑑x𝑑y\mu _{XY}(x,y)\mathrm{log}\frac{\mu _{XY}(x,y)}{\mu _X(x)\mu _Y(y)},$$
(6)
is non-negative and invariant under $`x\varphi (x)`$ and $`y\psi (y)`$. It is (the limit of) a true information,

$$I(X,Y)=\underset{\mathrm{\Delta }0}{lim}[H_{\mathrm{bin}}(X)+H_{\mathrm{bin}}(Y)H_{\mathrm{bin}}(X,Y)].$$
(7)
In applications, one usually has the data available in form of a statistical sample. To estimate $`I(X,Y)`$ one starts from $`N`$ bivariate measurements $`(x_i,y_i),i=1,\mathrm{}N`$ which are assumed to be iid (independent identically distributed) realizations. There exist numerous algorithms to estimate $`I(X,Y)`$ and entropies. We shall use in the following the MI estimators proposed recently in Ref. \[Kraskov et al., 2003\], and we refer to this paper for a review of alternative methods.

### 2.2 Algorithmic Information Theory

In contrast to Shannon theory where the basic objects are random variables and entropies are average informations, algorithmic information theory deals with individual symbol strings and with the actual information needed to specify them. To “specify” a sequence $`X`$ means here to give the necessary input to a universal computer $`U`$, such that $`U`$ prints $`X`$ on its output and stops. The analogon to entropy, called here usually the complexity $`K(X)`$ of $`X`$, is the minimal length of an input which leads to the output $`X`$, for fixed $`U`$. It depends on $`U`$, but it can be shown that this dependence is weak and can be neglected in the limit when $`K(X)`$ is large \[Li and Vitanyi, 1997\].

Let us denote the concatenation of two strings $`X`$ and $`Y`$ as $`XY`$. Its complexity is $`K(XY)`$. It is intuitively clear that $`K(XY)`$ should be larger than $`K(X)`$ but cannot be larger than the sum $`K(X)+K(Y)`$. Finally, one expects that $`K(X|Y)`$, defined as the minimal length of a program printing $`X`$ when $`Y`$ is furnished as auxiliary input, is related to $`K(XY)K(Y)`$. Indeed, one can show \[Li and Vitanyi, 1997\] (again within correction terms which become irrelevant asymptotically) that

$$0K(X|Y)K(XY)K(Y)K(X).$$
(8)
Notice the close similarity with Shannon entropy.

The algorithmic information in $`Y`$ about $`X`$ is finally defined as

$$I_{\mathrm{alg}}(X,Y)=K(X)K(X|Y)K(X)+K(Y)K(XY).$$
(9)
Within the same additive correction terms, one shows that it is symmetric, $`I_{\mathrm{alg}}(X,Y)=I_{\mathrm{alg}}(Y,X)`$, and can thus serve as an analogon to mutual information.

From the halting theorem it follows that $`K(X)`$ is in general not computable. But one can easily give upper bounds. Indeed, the length of any input which produces $`X`$ (e.g. by spelling it out verbatim) is an upper bound. Improved upper bounds are provided by any file compression algorithm such as gnuzip or UNIX “compress”. Good compression algorithms will give good approximations to $`K(X)`$, and algorithms whose performance does not depend on the input file length (in particular since they do not segment the file during compression) will be crucial for the following.

### 2.3 MI-Based Distance Measures

Mutual information itself is a similarity measure in the sense that small values imply large “distances” in a loose sense. But it would be useful to modify it such that the resulting quantity is a metric in the strict sense, i.e. satisfies the triangle inequality. Indeed, the first such metric is well known \[Cover and Thomas, 1991\]: The quantity

$$d(X,Y)=H(X|Y)+H(Y|X)=H(X,Y)I(X,Y)$$
(10)
satisfies the triangle inequality, in addition to being non-negative and symmetric and to satisfying $`d(X,X)=0`$. The proof proceeds by first showing that for any $`Z`$

$$H(X|Y)H(X,Z|Y)H(X|Z)+H(Z|Y).$$
(11)
But $`d(X,Y)`$ is not appropriate for our purposes. Since we want to compare the proximity between two single objects and that between two clusters containing maybe many objects, we would like the distance measure to be unbiased by the sizes of the clusters. As argued forcefully in \[Li et al., 2001, Li et al., 2002\], this is not true for $`I_{\mathrm{alg}}(X,Y)`$, and for the same reasons it is not true for $`I(X,Y)`$ or $`d(X,Y)`$ either: A mutual information of thousand bits should be considered as large, if $`X`$ and $`Y`$ themselves are just thousand bits long, but it should be considered as very small, if $`X`$ and $`Y`$ would each be huge, say one million bits.

As shown in \[Li et al., 2001, Li et al., 2002\] within the algorithmic framework, one can form two different distances which measure relative distance, i.e. which are normalized by dividing by a total entropy. We sketch here only the theorems and proofs for the Shannon version, they are indeed very similar to their algorithmic analoga in \[Li et al., 2001, Li et al., 2002\].

Theorem 1: The quantity

$$D(X,Y)=1\frac{I(X,Y)}{H(X,Y)}=\frac{d(X,Y)}{H(X,Y)}$$
(12)
is a metric, with $`D(X,X)=0`$ and $`D(X,Y)1`$ for all pairs $`(X,Y)`$.

Proof: Symmetry, positivity and boundedness are obvious. Since $`D(X,Y)`$ can be written as

$$D(X,Y)=\frac{H(X|Y)}{H(X,Y)}+\frac{H(Y|X)}{H(Y,X)},$$
(13)
it is sufficient for the proof of the triangle inequality to show that each of the two terms on the r.h.s. is bounded by an analogous inequality, i.e.

$$\frac{H(X|Y)}{H(X,Y)}\frac{H(X|Z)}{H(X,Z)}+\frac{H(Z|Y)}{H(Z,Y)}$$
(14)
and similarly for the second term. Eq.(14) is proven straightforwardly, using Eq.(11) and the basic inequalities $`H(X)0`$, $`H(X,Y)H(X,Y,Z)`$ and $`H(X|Z)0`$:

$`{\displaystyle \frac{H(X|Y)}{H(X,Y)}}`$ $`=`$ $`{\displaystyle \frac{H(X|Y)}{H(Y)+H(X|Y)}}`$ (15)
$``$ $`{\displaystyle \frac{H(X|Z)+H(Z|Y)}{H(Y)+H(X|Z)+H(Z|Y)}}`$
$`=`$ $`{\displaystyle \frac{H(X|Z)+H(Z|Y)}{H(X|Z)+H(Y,Z)}}`$
$``$ $`{\displaystyle \frac{H(X|Z)}{H(X|Z)+H(Z)}}+{\displaystyle \frac{H(Z|Y)}{H(Y,Z)}}`$
$`=`$ $`{\displaystyle \frac{H(X|Z)}{H(X,Z)}}+{\displaystyle \frac{H(Z|Y)}{H(Z,Y)}}.`$
Theorem 2: The quantity

$`D^{}(X,Y)`$ $`=`$ $`1{\displaystyle \frac{I(X,Y)}{\mathrm{max}\{H(X),H(Y)\}}}`$ (16)
$`=`$ $`{\displaystyle \frac{\mathrm{max}\{H(X|Y),H(Y|X)\}}{\mathrm{max}\{H(X),H(Y)\}}}`$
is also a metric, also with $`D^{}(X,X)=0`$ and $`D^{}(X,Y)1`$ for all pairs $`(X,Y)`$. It is sharper than $`D`$, i.e. $`D^{}(X,Y)D(X,Y)`$.

Proof: Again we have only to prove the triangle inequality, the other parts are trivial. For this we have to distinguish different cases \[Li et al., 2002\].
Case 1: $`\mathrm{max}\{H(Z),H(Y)\}H(X)`$. Using Eq.(11) we obtain

$`D^{}(X,Y)`$ $`=`$ $`{\displaystyle \frac{H(X|Y)}{H(X)}}{\displaystyle \frac{H(X|Z)}{H(X)}}+{\displaystyle \frac{H(Z|Y)}{H(Y)}}`$ (17)
$`=`$ $`D^{}(X,Z)+D^{}(Z,Y).`$
Case 2: $`\mathrm{max}\{H(Z),H(X)\}H(Y)`$. This is completely analogous.
Case 3: $`H(X)H(Y)<H(Z)`$. We now have to show that

$`D^{}(X,Y)`$ $`=`$ $`{\displaystyle \frac{H(Y|X)}{H(Y)}}{\displaystyle \frac{H(Y|Z)+H(Z|X)}{H(Y)}}`$ (18)
$`\stackrel{\mathrm{?}}{}`$ $`D^{}(X,Z)+D^{}(Z,Y)`$
$`=`$ $`{\displaystyle \frac{H(Z|X)}{H(Z)}}+{\displaystyle \frac{H(Z|Y)}{H(Z)}}.`$
Indeed, if the r.h.s. of the first line is less than 1, then

$`{\displaystyle \frac{H(Y|X)}{H(Y)}}`$ $``$ $`{\displaystyle \frac{H(Y|Z)+H(Z|X)}{H(Y)}}`$ (19)
$``$ $`{\displaystyle \frac{H(Y|Z)+H(Z|X)+H(Z)H(Y)}{H(Z)}}`$
$`=`$ $`{\displaystyle \frac{H(Z|Y)+H(Z|X)}{H(Z)}},`$
and Eq.(18) holds. If it is larger than 1, then also $`(H(Z|Y)+H(Z|X))/H(Z)1`$. Eq.(18) must now also hold, since $`H(Y|X)/H(Y)1`$.
Case 4: $`H(Y)H(X)<H(Z)`$. This is completely analogous to case 3.

Apart from scaling correctly with the total information, in contrast to $`d(X,Y)`$, the algorithmic analog to $`D(X,Y)`$ and $`D^{}(X,Y)`$ are also universal \[Li et al., 2002\]. Essentially this means that if $`XY`$ according to any non-trivial distance measure, then $`XY`$ also according to $`D`$, and even more so (by factor up to 2) according to $`D^{}`$. In contrast to the other properties of $`D`$ and $`D^{}`$, this is not easy to carry over from algorithmic to Shannon theory. The proof in Ref. \[Li et al., 2002\] depends on $`X`$ and $`Y`$ being discrete, which is obviously not true for probability distributions. Based on the universality argument, it was argued in \[Li et al., 2002\] that $`D^{}`$ should be superior to $`D`$, but the numerical studies shown in that reference did not show a clear difference between them. In the following we shall therefore use primarily $`D`$ for simplicity, but we checked that using $`D^{}`$ did not give systematically better results.

A major difficulty appears in the Shannon framework, if we deal with continuous random variables. As we mentioned above, Shannon informations are only finite for coarse-grained variables, while they diverge if the resolution tends to zero. This means that dividing MI by the entropy as in the definitions of $`D`$ and $`D^{}`$ becomes problematic. One has essentially two alternative possibilities. The first is to actually introduce some coarse-graining, although it would have been necessary for the definition of $`I(X,Y)`$, and divide by the coarse-grained entropies. This introduces an arbitrariness, since the scale $`\mathrm{\Delta }`$ is completely ad hoc, unless it can be fixed by some independent arguments. We have found no such arguments, and thus we propose the second alternative. There we take $`\mathrm{\Delta }0`$. In this case $`H(X)m_x\mathrm{log}\mathrm{\Delta }`$, with $`m_x`$ being the dimension of $`X`$. In this limit $`D`$ and $`D^{}`$ would tend to 1. But using similarity measures

$`S(X,Y)=(1D(X,Y))\mathrm{log}(1/\mathrm{\Delta }),`$ (20)
$`S^{}(X,Y)=(1D^{}(X,Y))\mathrm{log}(1/\mathrm{\Delta })`$ (21)
instead of $`D`$ and $`D^{}`$ gives exactly the same results in MIC, and

$$S(X,Y)=\frac{I(X,Y)}{m_x+m_y},S^{}(X,Y)=\frac{I(X,Y)}{\mathrm{max}\{m_x,m_y\}}.$$
(22)
Thus, when dealing with continuous variables, we divide the MI either by the sum or by the maximum of the dimensions. When starting with scalar variables and when $`X`$ is a cluster variable obtained by joining $`m`$ elementary variables, then its dimension is just $`m_x=m`$.

## 3 A Phylogenetic Tree for Mammals

As a first application, we study the mitochondrial DNA of a group of 34 mammals (see Fig. 1). Exactly the same data \[Gen, a\] had previously been analyzed in \[Li et al., 2001, Reyes et al., 2000\]. This group includes among others222opossum (Didelphis virginiana), wallaroo (Macropus robustus), and platypus (Ornithorhyncus anatinus) some rodents333rabbit (Oryctolagus cuniculus), guinea pig (Cavia porcellus), fat dormouse (Glis glis), rat (Rattus norvegicus), squirrel (Scuirus vulgaris), and mouse (Mus musculus), ferungulates444horse (Equu caballus), donkey (Equus asinus), Indian rhinoceros (Rhinoceros unicornis), white rhinoceros (Ceratotherium simum), harbor seal (Phoca vitulina), grey seal (Halichoerus grypus), cat (Felis catus), dog (Canis familiaris), fin whale (Balenoptera physalus), blue whale (Balenoptera musculus), cow (Bos taurus), sheep (Ovis aries), pig (Sus scrofa), hippopotamus (Hippopotamus amphibius), neotropical fruit bat (Artibeus jamaicensis), African elephant (Loxodonta africana), aardvark (Orycteropus afer), and armadillo (Dasypus novemcintus), and primates555human (Homo sapiens), common chimpanzee (Pan troglodytes), pigmy chimpanzee (Pan paniscus), gorilla (Gorilla gorilla), orangutan (Pongo pygmaeus), gibbon (Hylobates lar), and baboon (Papio hamadryas). It had been chosen in \[Li et al., 2001\] because of doubts about the relative closeness among these three groups \[Cao et al., 1998, Reyes et al., 2000\].

Obviously, we are here dealing with the algorithmic version of information theory, and informations are estimated by lossless data compression. For constructing the proximity matrix between individual taxa, we proceed essentially a in Ref. \[Li et al., 2001\]. But in addition to using the special compression program GenCompress \[Gen, b\], we also tested several general purpose compression programs such as BWTzip \[BWT,\] and the UNIX tool bzip2.

In Ref. \[Li et al., 2001\], this proximity matrix was then used as the input to a standard HC algorithm (neighbour-joining and hypercleaning) to produce an evolutionary tree. It is here where our treatment deviates crucially. We used the MIC algorithm described in Sec. 1, with distance $`D(X,Y)`$. The joining of two clusters (the third step in the MIC algorithm) is obtained by simply concatenating the DNA sequences. There is of course an arbitrariness in the order of concatenation sequences: $`XY`$ and $`YX`$ give in general compressed sequences of different lengths. But we found this to have negligible effect on the evolutionary tree. The resulting evolutionary tree obtained with Gencompress is shown in Fig. 1.

As shown in Fig. 1 the overall structure of this tree closely resembles the one shown in Ref. \[Reyes et al., 2000\]. All primates are correctly clustered and also the relative order of the ferungulates is in accordance with Ref. \[Reyes et al., 2000\]. On the other hand, there are a number of connections which obviously do not reflect the true evolutionary tree, see for example the guinea pig with bat and elephant with platypus. But the latter two, inspite of being joined together, have a very large distance from each other, thus their clustering just reflects the fact that neither the platypus nor the elephant have other close relatives in the sample. All in all, however, already the results shown in Fig. 1 capture surprisingly well the overall structure shown in Ref. \[Reyes et al., 2000\]. Dividing MI by the total information is essential for this success. If we had used the non-normalized $`I_{\mathrm{alg}}(X,Y)`$ itself, the clustering algorithm used in \[Li et al., 2001\] would not change much, since all 34 DNA sequences have roughly the same length. But our MIC algorithm would be completely screwed up: After the first cluster formation, we have DNA sequences of very different lengths, and longer sequences tend also to have larger MI, even if they are not closely related.

A heuristic reasoning for the use of MIC for the reconstruction of an evolutionary tree might be given as follows: Suppose that a proximity matrix has been calculated for a set of DNA sequences and the smallest distance is found for the pair $`(X,Y)`$. Ideally, one would remove the sequences $`X`$ and $`Y`$, replace them by the sequence of the common ancestor (say $`Z`$) of the two species, update the proximity matrix to find the smallest entry in the reduced set of species, and so on. But the DNA sequence of the common ancestor is not available. One solution might be that one tries to reconstruct it by making some compromise between the sequences $`X`$ and $`Y`$. Instead, we essentially propose to concatenate the sequences $`X`$ and $`Y`$. This will of course not lead to a plausible sequence of the common ancestor, but it will optimally represent the information about the common ancestor. During the evolution since the time of the ancestor $`Z`$, some parts of its genome might have changed both in $`X`$ and in $`Y`$. These parts are of little use in constructing any phylogenetic tree. Other parts might not have changed in either. They are recognized anyhow by any sensible algorithm. Finally, some parts of its genome will have mutated significantly in $`X`$ but not in $`Y`$, and vice versa. This information is essential to find the correct way through higher hierarchy levels of the evolutionary tree, and it is preserved in concatenating.

## 4 Clustering of Minimally Dependent Components in an Electrocardiogram

As our second application we choose a case where Shannon theory is the proper setting. We show in Fig. 2 an ECG recorded from the abdomen and thorax of a pregnant woman \[De Moor, 1997\] (8 channels, sampling rate 500 Hz, 5 s total). It is already seen from this graph that there are at least two important components in this ECG: the heartbeat of the mother, with a frequency of $`3`$ beat/s, and the heartbeat of the fetus with roughly twice this frequency. Both are not synchronized. In addition there is noise from various sources (muscle activity, measurement noise, etc.). While it is easy to detect anomalies in the mother’s ECG from such a recording, it would be difficult to detect them in the fetal ECG.

As a first approximation we can assume that the total ECG is a linear superposition of several independent sources (mother, child, noise1, noise2,…). A standard method to disentangle such superpositions is independent component analysis (ICA) \[Hyvärinen et al., 2001\]. In the simplest case one has $`n`$ independent sources $`s_i(t),i=1\mathrm{}n`$ and $`n`$ measured channels $`x_i(t)`$ obtained by instantaneous superpositions with a time independent non-singular matrix $`𝐀`$,

$$x_i(t)=\underset{j=1}{\overset{n}{}}A_{ij}s_j(t).$$
(23)
In this case the sources can be reconstructed by applying the inverse transformation $`𝐖=𝐀^1`$ which is obtained by minimizing the (estimated) mutual informations between the transformed components $`y_i(t)=_{j=1}^nW_{ij}x_j(t)`$. If some of the sources are Gaussian, this leads to ambiguities \[Hyvärinen et al., 2001\], but it gives a unique solution if the sources have more structure.

In reality things are not so simple. For instance, the sources might not be independent, the number of sources (including noise sources!) might be different from the number of channels, and the mixing might involve delays. For the present case this implies that the heartbeat of the mother is seen in several reconstructed components $`y_i`$, and that the supposedly “independent” components are not independent at all. In particular, all components $`y_i`$ which have large contributions from the mother form a cluster with large intra-cluster MIs and small inter-cluster MIs. The same is true for the fetal ECG, albeit less pronounced. It is thus our aim to
1) optimally decompose the signals into least dependent components;
2) cluster these components hierarchically such that the most dependent ones are grouped together;
3) decide on an optimal level of the hierarchy, such that the clusters make most sense physiologically;
4) project onto these clusters and apply the inverse transformations to obtain cleaned signals for the sources of interest.

Technically we proceeded as follows \[Stögbauer et al., 2004\]:

Since we expect different delays in the different channels, we first used Takens delay embedding \[Takens, 1980\] with time delay 0.002 s and embedding dimension 3, resulting in $`24`$ channels. We then formed 24 linear combinations $`y_i(t)`$ and determined the de-mixing coefficients $`W_{ij}`$ by minimizing the overall mutual information between them, using the MI estimator proposed in \[Kraskov et al., 2003\]. There, two classes of estimators were introduced, one with square and the other with rectangular neighbourhoods. Within each class, one can use the number of neighbours, called $`k`$ in the following, on which the estimate is based. Small values of $`k`$ lead to a small bias but to large statistical errors, while the opposite is true for large $`k`$. But even for very large $`k`$ the bias is zero when the true MI is zero, and it is systematically such that absolute values of the MI are underestimated. Therefore this bias does not affect the determination of the optimal de-mixing matrix. But it depends on the dimension of the random variables, therefore large values of $`k`$ are not suitable for the clustering. We thus proceeded as follows: We first used $`k=100`$ and square neighbourhoods to obtain the least dependent components $`y_i(t)`$, and then used $`k=3`$ with rectangular neighbourhoods for the clustering. The resulting least dependent components are shown in Fig. 3. They are sorted such that the first components (1 - 5) are dominated by the mother’s ECG, while the next three contain large contributions from the fetus. The rest contains mostly noise, although some seem to be still mixed.

These results obtained by visual inspection are fully supported by the cluster analysis. The dendrogram is shown in Fig. 4. In constructing it we used $`S(X,Y)`$ (Eq.(22)) as similarity measure to find the correct topology. Again we would have obtained much worse results if we had not normalized it by dividing MI by $`m_X+m_Y`$. In plotting the actual dendrogram, however, we used the MI of the cluster to determine the height at which the two daughters join. The MI of the first five channels, e.g., is $`1.43`$, while that of channels 6 to 8 is $`0.34`$. For any two clusters (tuples) $`X=X_1\mathrm{}X_n`$ and $`Y=Y_1\mathrm{}Y_m`$ one has $`I(X,Y)I(X)+I(Y)`$. This guarantees, if the MI is estimated correctly, that the tree is drawn properly. The two slight glitches (when clusters (1–14) and (15–18) join, and when (21–22) is joined with 23) result from small errors in estimating MI. They do in no way effect our conclusions.

In Fig. 4 one can clearly see two big clusters corresponding to the mother and to the child. There are also some small clusters which should be considered as noise. For reconstructing the mother and child contributions to Fig. 2, we have to decide on one specific clustering from the entire hierarchy. We decided to make the cut such that mother and child are separated. The resulting clusters are indicated in Fig. 4 and were already anticipated in sorting the channels. Reconstructing the original ECG from the child components only, we obtain Fig. 5.

## 5 Conclusions

We have shown that MI can not only be used as a proximity measure in clustering, but that it also suggests a conceptually very simple and natural hierarchical clustering algorithm. We do not claim that this algorithm, called mutual information clustering (MIC), is always superior to other algorithms. Indeed, MI is in general not easy to estimate. Obviously, when only crude estimates are possible, also MIC will not give optimal results. But as MI estimates are becoming better, also the results of MIC should improve. The present paper was partly triggered by the development of a new class of MI estimators for continuous random variables which have very small bias and also rather small variances \[Kraskov et al., 2003\].

We have illustrated our method with two applications, one from genetics and one from cardiology. For neither application MIC might give the very best clustering, but it seems promising that one common method gives decent results for both, although they are very different.

The results of MIC should improve, if more data become available. This is trivial, if we mean by that longer time sequences in the application to ECG, and longer parts of the genome in the application of Sec.3. It is less trivial that we expect MIC to make fewer mistakes in a phylogenetic tree, when more species are included. The reason is that close-by species will be correctly joined anyhow, and families – which now are represented only by single species and thus are poorly characterized – will be much better described by the concatenated genomes if more species are included.

There are two versions of information theory, algorithmic and probabilistic, and therefore there are also two variants of MI and of MIC. We discussed in detail one application of each, and showed that indeed common concepts were involved in both. In particular it was crucial to normalize MI properly, so that it is essentially the relative MI which is used as proximity measure. For conventional clustering algorithms using algorithmic MI as proximity measure this had already been stressed in \[Li et al., 2001, Li et al., 2002\], but it is even more important in MIC, both in the algorithmic and in the probabilistic versions.

In the probabilistic version, one studies the clustering of probability distributions. But usually distributions are not provided as such, but are given implicitly by finite random samples drawn (more or less) independently from them. On the other hand, the full power of algorithmic information theory is only reached for infinitely long sequences, and in this limit any individual sequence defines a sequence of probability measures on finite subsequences. Thus the strict distinction between the two theories is somewhat blurred in practice. Nevertheless, one should not confuse the similarity between two sequences (two English books, say) and that between their subsequence statistics. Two sequences are maximally different if they are completely random, but their statistics for short subsequences is then identical (all subsequences appear in both with equal probabilities). Thus one should always be aware of what similarities or independencies one is looking for. The fact that MI can be used in similar ways for all these problems is not trivial.

We would like to thank Arndt von Haesseler, Walter Nadler and Volker Roth for many useful discussions.


text:
# Currents Induced by Charges Moving in Semiconductor

## Introduction

The method of computation of currents induced on electrodes by charges moving in vacuum was introduced in , . In presented paper, this method is generalised for charges moving in a media.

## Derivation of the Induced Current

Consider the moving charge, $`q_m`$, inside a semiconductor and any number of conductors at constant potentials for one of which, say A (following Ramo’s approach ), the induced current is desired. The potential of the electrostatic field inside the semiconductor under static or low frequency conditions  satisfies

$$\epsilon \epsilon _0^2V=\rho ,$$
(1)
where $`\rho `$ is a charge density, $`\epsilon _0`$ is permittivity of free space and $`\epsilon `$ is permittivity of the media.

The boundaries are surfaces of electrodes and tiny equipotential sphere surrounding moving charge. The integral over a boundary surface $`S_i`$ equals to

$$_{S_i}\epsilon \epsilon _0\frac{V}{n}𝑑s=Q_i,$$
(2)
where $`n`$ is the outward normal to the surface and $`Q_i`$ is the charge on the i-th electrode or moving charge.

Also consider the same geometry of conductors in the free space but without moving charge. Potential of the electrode A is rised to unit potential and all other conductors are grounded. Potential in this case, $`V^{}`$, satisfies $`^2V^{}=0`$ everywhere between conductors.

Applying Green’s theorem to the vector function

$$𝐟=\epsilon \epsilon _0(V^{}VVV^{}),$$
(3)
we obtain the following equation

$$_{Volume}\epsilon \epsilon _0(V^{}^2VV^2V^{})𝑑v=_{Surface}\epsilon \epsilon _0(V^{}\frac{V}{n}V\frac{V^{}}{n})𝑑s$$
(4)
Using, that by design, $`V^{}=0`$ on the surfaces of all electrodes but A and $`^2V^{}=0`$ in the space between boundaries Eq. 4 can be rewritten as

$$_{Volume}\epsilon \epsilon _0V^{}^2Vdv=Q_A+V_m^{}q_m\underset{electrodes}{}\epsilon V_iQ_i^{}$$
(5)
or

$$Q_A=q_mV_m^{}+\underset{electrodes}{}\epsilon V_iQ_i^{}+_{Volume}V^{}\epsilon \epsilon _0^2Vdv$$
(6)
$$i_A=\frac{dQ_a}{dt}=q_m\frac{dV_m^{}}{dt}+\frac{}{t}_{Volume}V^{}\epsilon \epsilon _0^2Vdv$$
(7)
This can be revritten as

$$i_A=q_m\frac{dV_m^{}}{d𝐫}\frac{d𝐫}{dt}+_{Volume}V^{}\frac{\rho }{t}𝑑v$$
(8)
where $`𝐫`$ is the direction of the motion. In terms of weighting field, $`𝐄^{}=V^{}`$ this became

$$i_A=q_m𝐯𝐄^{}+_{Volume}V^{}\frac{\rho }{t}𝑑v$$
(9)
where the integration is over the volume of the media. In the case when media is vacuum or isolator, the second term is zero and Eq. 9 became Ramo-Shockley’s equation. In the case of semiconductor, this term is a weighted sum of currents generated by moving charge inside the semiconductor.