text:
# Comment on “Statistics of the Lyapunov Exponent in 1D Random periodic-on-Average Systems” [Phys. Rev. Lett. 81, 5390, 1998]

\[

\]

In this comment, we would like to point out possible numerical errors in a recent Letter by Deych et al..

In , the authors pointed out an interesting observation that in some 1D random periodic-on-average systems, states from pass bands and band edges of the underlying band structure follow single parameter scaling with universal behavior, whereas states from the interior of the band gaps do not have universal behavior and require two parameters to describe their scaling properties. When the degree of disorder exceeds a certain threshold the single parameter scaling is restored for an entire band gap. Upon careful inspection, we find that while these conclusions are qualitatively valid, there may be some numerical errors in the results of .

The model and all the parameters used in the following computation are the same as in . The model in  is essentially the same as the second model in . First we compute the band structure using Eq. (6) of  the result is presented in Fig. 1. It shows that the first band gap occurs between $`k=1.456`$ and $`1.543`$ and there is a very narrow second gap located between $`k=2.992`$ and $`3.005`$. The third band gap is at around $`k=4.5`$. In , however, the results indicates a band gap at around $`k=2.86`$ which is absent in our computation. To validate our program, we used the parameters from  in the computation and obtained the same band structure results as in . Moreover, by adjusting the parameters we were also able to reproduce almost exactly the same result for the disorder dependence of the Lyapunov exponent as in Fig. 1 of .

Following , we then compute the disorder dependence of the Lyapunov exponent on the standard deviation in the randomness of the layers’ thickness denoted by $`\sigma `$, and the frequency dependence of the Lyapunov exponent and its variance for frequencies covering the first three band gaps. The results are shown in Fig. 2. It is clear that the results for frequencies covering the first band gaps are qualitatively in agreement with the results in Fig. 1 and Fig. 2 of , except that the band gap in  is implied at around $`k=2.86`$. In particular, as the randomness increases, all the Lyapunov exponent curves merge. For frequencies in the vicinity of the second band gap, however, the exponent curves come across into each other at a certain level of disorder, but split again as the amount of disorder increases further. This is shown in Fig. 2(a2). In this case, for the same level of the disorder the double peaked structure of the variance observed in the first band gap disappears (Fig. 1(b2)). When the disorder decreases, to $`\sigma =0.00577`$ say, the double peaks appear again. For frequencies near the third band gap, as indicated by Fig. 2(a3) and (b3) the results appear to follow the features at frequencies around the first band gap. For comparison, we also plot in Fig. 3 the dependence of the Lyapunov exponent on the standard deviation $`\sigma `$ of the layer randomness at the frequencies indicated in Fig. 2 of . At these frequencies, all curves seem to merge.

This work received support from National Science Council (Grant No. NSC-89-2112-M008-008 and NSC-89-2611-M008-002).


format: text

text:
# References

In Ref. , we argued that dilatonic domain walls with $`\mathrm{\Delta }<2`$ (see just below Eq. (5) for definition of $`\mathrm{\Delta }`$) localize gravity, since the potential term in the Schrödinger-type equation satisfied by the metric fluctuation is volcano-shaped just like that of the Randall-Sundrum (RS) domain wall . In the subsequent paper , we found out that gravity can be localized around the positive tension domain wall with $`\mathrm{\Delta }<0`$, since the effective gravitational constant in one lower dimensions is nonzero for such case. Then, later it is realized  that the graviton can be localized on the dilatonic domain wall for any values of $`\mathrm{\Delta }`$ as long as the tension of the domain wall is positive, since for such case the normalizable Kaluza-Klein (KK) zero mode exists and the effective gravitational constant in one lower dimensions is nonzero. We feel that these different (but expanded) results scattered around our previous works may cause confusion among readers. Furthermore, perhaps to worsen the confusion, in Ref. , which later reproduces these our previous papers, it is argued that gravity can be trapped on the dilatonic domain walls only for the $`\mathrm{\Delta }2`$ case, since only for such case normalizable graviton KK mode exists. It is a purpose of this brief note to gather the results on the localization of gravity on the dilatonic domain wall, which are scattered around our different previous papers, to give unified and coherent discussion in order to prevent misleading readers, and hopefully to clarify the contradicting result of Ref. . We also elaborate on the possible embeddings of the RS type brane world scenario into string theories through dilatonic and non-dilatonic domain walls, which follows straightforwardly from our previous works but was not elaborated explicitly.

The action for the dilatonic domain wall with an arbitrary dilaton coupling parameter $`a`$ in arbitrary spacetime dimensions $`D`$ is

$$S_{\mathrm{bulk}}=\frac{1}{2\kappa _D^2}d^Dx\sqrt{G}\left[_G\frac{4}{D2}_M\varphi ^M\varphi +e^{2a\varphi }\mathrm{\Lambda }\right].$$
(1)
The worldvolume action for the dilatonic domain wall with energy density (or tension) $`\sigma _{DW}`$ is

$$S_{\mathrm{DW}}=d^{D1}x\sqrt{\gamma }\sigma _{DW}e^{a\varphi },$$
(2)
where $`\gamma `$ is the determinant of the induced metric $`\gamma _{\mu \nu }=_\mu X^M_\nu X^NG_{MN}`$ on the domain wall worldvolume.

The dilatonic domain wall solution in the standard warp factor form is :

$`G_{MN}dx^Mdx^N=𝒲\left[dt^2+dx_1^2+\mathrm{}+dx_{D2}^2\right]+dy^2,`$ (3)
$`\varphi ={\displaystyle \frac{1}{a}}\mathrm{ln}(1+Ky),𝒲=(1+Ky)^{\frac{8}{(D2)^2a^2}},`$ (4)
where the constant $`K`$ can take following positive or negative value independently on both sides of the domain wall:

$$K=\pm \frac{(D2)a^2}{2}\sqrt{\frac{D2}{4(D1)a^2(D2)^2}\mathrm{\Lambda }}=\pm \frac{(D2)a^2}{2}\sqrt{\frac{\mathrm{\Lambda }}{2\mathrm{\Delta }}},$$
(5)
where $`\mathrm{\Delta }(D2)a^2/22(D1)/(D2)`$. From this expression for $`K`$, one can see that the domain solution of the form (4) exists only for $`\mathrm{\Delta }<0`$ \[$`\mathrm{\Delta }>0`$\] when $`\mathrm{\Lambda }>0`$ \[$`\mathrm{\Lambda }<0`$\]. The boundary condition at $`y=0`$ fixes the domain wall tension $`\sigma _{DW}`$ to be the following fine-tuned form :

$$\sigma _{DW}=\frac{1}{\kappa _D^2}\frac{4}{(D2)a^2}(K_{}K_+),$$
(6)
where $`K_{}`$ \[$`K_+`$\] denotes the value of $`K`$ at $`y<0`$ \[$`y>0`$\]. So, $`\sigma _{DW}`$ is positive \[negative\] when $`K_{}>0`$ and $`K_+<0`$ \[$`K_{}<0`$ and $`K_+>0`$\], in which case there are naked singularities at finite nonzero $`y`$ \[no singularity at $`y0`$\]. And $`\sigma _{DW}`$ is zero when $`K_+`$ and $`K_{}`$ have the same sign. The $`(D1)`$-dimensional gravitational constant has non-zero value given by $`\kappa _{D1}^2=\frac{\mathrm{\Delta }+4}{2}\sqrt{\frac{\mathrm{\Lambda }}{2\mathrm{\Delta }}}\kappa _D^2`$ when $`\sigma _{DW}>0`$; otherwise, $`\kappa _{D1}^2=0`$ . So, a necessary condition for localizing gravity on the domain wall is $`\sigma _{DW}>0`$, which is possible for any values of $`\mathrm{\Delta }`$.

By redefining the transverse coordinate of the domain wall, one can put the solution into the following conformally flat form :

$`G_{MN}dx^Mdx^N=𝒞\left[dt^2+dx_1^2+\mathrm{}+dx_{D2}^2+dz^2\right],`$ (7)
$`\varphi ={\displaystyle \frac{(D2)a}{2(\mathrm{\Delta }+2)}}\mathrm{ln}(1+\overline{K}z),𝒞=(1+\overline{K}z)^{\frac{4}{(D2)(\mathrm{\Delta }+2)}},`$ (8)
where the constant $`\overline{K}`$ is given by

$$\overline{K}=\eta \frac{(D2)^2a^24}{(D2)^2a^2}K=\pm \eta (\mathrm{\Delta }+2)\sqrt{\frac{\mathrm{\Lambda }}{2\mathrm{\Delta }}},$$
(9)
where $`\pm `$ is the same as $`\pm `$ in Eq. (5) and the sign ambiguity $`\eta =\pm `$ resulting from the coordinate transformation to the conformally flat metric can be fixed to be $`\eta =+`$ by demanding that $`\kappa _{D1}^2`$ should be calculated to take the same values as above. From now on, we consider only the $`𝐙_2`$-symmetric solution with $`\sigma _{DW}>0`$, for which $`\kappa _{D1}0`$. So, in the remainder of this paper, we let $`𝒲=(1+K|y|)^{\frac{8}{(D2)^2a^2}}`$ and $`𝒞=(1+\overline{K}|z|)^{\frac{4}{(D2)(\mathrm{\Delta }+2)}}`$ with $`K`$ and $`\overline{K}`$ respectively given by Eqs. (5) and (9) with $``$ sign. \[For the choice of $`+`$ sign, $`\sigma _{DW}<0`$ and $`\kappa _{D1}=0`$.\]

One can determine whether gravity can be localized on the domain wall also by studying the linearized Einstein equations satisfied by the metric fluctuation $`h_{\mu \nu }(x^\rho ,z)=\widehat{h}_{\mu \nu }^{(m)}(x^\rho )𝒞^{(D2)/4}\psi _m(z)`$ in the RS gauge. When $`\eta ^{\mu \nu }_\mu _\nu \widehat{h}_{\mu \nu }^{(m)}=m^2\widehat{h}_{\mu \nu }^{(m)}`$, $`\psi _m(z)`$ satisfies the following Schrödinger-type equation :

$$\frac{d^2\psi _m}{dz^2}+V(z)\psi _m=m^2\psi _m,$$
(10)
with the potential given by

$$V(z)=\frac{D2}{16}\left[(D6)\left(\frac{𝒞^{}}{𝒞}\right)^2+4\frac{𝒞^{\prime \prime }}{𝒞}\right]=\frac{(1+\mathrm{\Delta })\mathrm{\Lambda }}{2\mathrm{\Delta }(1+\overline{K}|z|)^2}+\frac{2\overline{K}}{\mathrm{\Delta }+2}\delta (z).$$
(11)
Note, the potential expression in Eq. (24) of Ref.  is related to this expression through $`Q=\frac{\mathrm{\Delta }}{\mathrm{\Delta }+2}\overline{K}=\pm \mathrm{\Delta }\sqrt{\frac{\mathrm{\Lambda }}{2\mathrm{\Delta }}}`$ with the sign $`\pm `$ being that in Eq. (9). \[In Refs. , we considered $`Q>0`$ case, only, namely $``$ sign choice ($`\sigma _{DW}>0`$) for $`\mathrm{\Delta }<0`$ and $`+`$ sign choice ($`\sigma _{DW}<0`$) for $`\mathrm{\Delta }>0`$. In this paper, we consider more general case, as we did in Ref. .\] Note, $`V(z)`$ has attractive \[repulsive\] $`\delta `$-function potential term when $`\sigma _{DW}>0`$ \[$`\sigma _{DW}<0`$\] for any value of $`\mathrm{\Delta }`$. So, for any value of $`\mathrm{\Delta }`$, the graviton KK zero mode $`\psi _0𝒞^{(D2)/4}`$, satisfying the boundary condition $`\psi _0^{}(0^+)\psi _0^{}(0^{})=\frac{2\overline{K}}{\mathrm{\Delta }+2}\psi _0(0)`$, is normalizable (i.e., $`𝑑z|\psi _0|^2<\mathrm{}`$) 222Note, for the positive tension domain walls, $`\overline{K}`$ is positive \[negative\] for $`\mathrm{\Delta }<2`$ \[$`\mathrm{\Delta }>2`$\] and therefore the metric (8) is well-defined over the interval $`\mathrm{}<z<\mathrm{}`$ \[$`\overline{K}^1z\overline{K}^1`$\], over which the normalization integral is integrated, as long as we choose $``$ sign in Eq. (9) . Even for the $`\sigma _{DW}<0`$ case, the graviton KK zero mode is still given by $`\psi _0𝒞^{(D2)/4}`$ for any values of $`\mathrm{\Delta }`$ but is not normalizable. Particularly, for the $`\mathrm{\Delta }<2`$ case with $`\sigma _{DW}>0`$, $`V(z)`$ resembles the volcano potential of the RS model  and therefore has the same KK mode spectrum structure and the same type of correction to the Newtonian gravitational potential from the massive KK modes as the RS domain wall case, as was noted in Ref. . However, note that even for the $`\mathrm{\Delta }>2`$ case with $`\sigma _{DW}>0`$, the normalizable KK zero mode exists, even if $`V(z)`$ does not take the volcano potential form. The structure of the massive KK mode spectrum for the $`\mathrm{\Delta }>2`$ case is discussed in Ref. . \[Note, the possible tachyonic modes for the $`\mathrm{\Delta }>1`$ case mentioned in Ref.  do not exist, since the result of Ref.  ensures that $`m^20`$. For the $`\mathrm{\Delta }>0`$ case, Ref.  considers the case of repulsive $`\delta `$-function potential (i.e., $`+`$ sign choice in Eq. (9)), but the structure of the massive KK mode spectrum discussed in Ref.  still remains the same even for the attractive $`\delta `$-function potential case.\] However, in the $`\mathrm{\Delta }>0`$ case (i.e., the $`\mathrm{\Lambda }<0`$ case), the $`(D1)`$-dimensional effective action has diverging cosmological constant term . It is found out  that to avoid such divergence one has to cut off the transverse space by introducing additional domain wall with the fine-tuned tension between $`z=0`$ and the curvature singularity.

To sum up, the normalizable graviton KK zero mode exists and the effective gravitational constant in one lower dimension is nonzero for any values of $`\mathrm{\Delta }`$ when $`\sigma _{DW}>0`$, but the introduction of additional domain wall is required for the $`\mathrm{\Delta }>0`$ case to remove diverging cosmological constant term in the effective action in one lower dimensions. And when $`\sigma _{DW}0`$, the graviton KK zero mode is not normalizable and the gravitational constant in one lower dimension is zero for any values of $`\mathrm{\Delta }`$. \[Note, the graviton KK modes for the $`\mathrm{\Delta }=2`$ case was not studied by us, but was later studied in Ref. . It is shown there that the normalizable graviton KK zero mode exists and there is the continuum of the massive KK mode with mass gap for the $`\mathrm{\Delta }=2`$ case with $`\sigma _{DW}>0`$.\]

We now comment on the results of Ref.  which contradict the results of our previous papers. The authors of Ref.  claim that the normalizable graviton KK zero mode exists only for the $`\mathrm{\Delta }2`$ case with $`\sigma _{DW}>0`$, for which $`\overline{K}`$ (which corresponds to $`k`$ in Ref. ) is positive, and the introduction of additional domain wall is necessary for the $`\mathrm{\Delta }>2`$ case in order to trap gravity. However, our careful analysis shows that gravity can be trapped even for the $`\mathrm{\Delta }>2`$ case, for which $`k<0`$ for the positive tension domain wall, and the introduction of additional domain wall is necessary for the $`\mathrm{\Delta }>0`$ case not to trap gravity but to remove diverging cosmological constant term in the $`(D1)`$-dimensional effective action. As necessary conditions for trapping gravity, Ref.  states that the conformal factor $`𝒞`$ should vanish at large $`|z|`$ and the $`\delta `$-function source should have a positive tension so that the potential term in the Schrödinger equation can be volcano-like. However, Ref.  fails to notice that $`𝒞`$ vanishes at finite $`|z|`$ (instead of at $`|z|=\mathrm{}`$) for the $`\mathrm{\Delta }>2`$ and $`\sigma _{DW}>0`$ case. And as we discussed, the potential term needs not be volcano-like in order to support the normalizable graviton KK zero mode.

We now discuss string theory embeddings of the above-discussed domain wall solutions which localize gravity. For this purpose, we begin by studying the uplifting of such domain walls to the dilatonic $`p`$-branes in $`D^{}>D`$ dimensions. The action for the dilatonic $`p`$-brane is given by

$$S_p=\frac{1}{2\kappa _D^{}^2}d^D^{}x\sqrt{\widehat{G}}\left[_{\widehat{G}}\frac{4}{D^{}2}(\phi )^2\frac{1}{2(p+2)!}e^{2a_p\phi }F_{p+2}^2\right].$$
(12)
The dilatonic $`p`$-brane solution is characterized by the parameter $`\mathrm{\Delta }_p=(D^{}2)a_p^2/2+2(p+1)(D^{}p3)/(D^{}2)`$. (See Eq. (2) of Ref.  for the explicit solution with this convention for the action.) Any single-charged branes in string theories are special cases of the dilatonic $`p`$-branes. The dilatonic $`p`$-branes can also be realized from any intersecting branes (with $`N`$ numbers of constituents) in string theories by setting the charges of the constituent branes equal to one another and then compactifying along the relative transverse directions (and possibly overall transverse and longitudinal directions). When the dilatonic $`p`$-branes are embedded in string theories in such a manner, the parameter $`\mathrm{\Delta }_p`$ takes only the special values $`4/N`$ and is invariant under the compactification on a Ricci flat manifold involving consistent truncation . So, the dilatonic domain wall solutions obtained from (intersecting) branes in string theories through the Scherk-Schwarz dimensional reduction  on a Ricci flat manifold have $`\mathrm{\Delta }=4/N`$, only. However, one can have dilatonic domain wall solutions in string theories with different values of $`\mathrm{\Delta }`$ through the compactification on spheres. By compactifying the $`D^{}`$-dimensional dilatonic $`p`$-brane (in the near horizon limit) with $`\mathrm{\Delta }_p=4/N`$ on $`S^{D^{}p2}`$, one obtains the dilatonic domain wall in $`(p+2)`$-dimensions with the dilaton coupling parameter $`a`$ and $`\mathrm{\Delta }`$ given by

$`|a|`$ $`=`$ $`{\displaystyle \frac{2}{p}}\sqrt{{\displaystyle \frac{2(D^{}2)N(p+1)(D^{}p3)}{2(D^{}p2)N(D^{}p3)}}},`$ (13)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{4(D^{}p3)}{2(D^{}p2)N(D^{}p3)}}.`$ (14)
By applying these results, we now elaborate on various possible string theory embeddings of dilatonic domain walls in five dimensions.

First, we discuss the case of the dilatonic domain walls with $`\mathrm{\Delta }=4/N`$ ($`N𝐙^+`$). To obtain such domain walls, we $`(i)`$ start with (intersecting) branes in ten or eleven dimensions (with equal constituent brane charges) with at least 3-dimensional (overall) longitudinal space and at least 1-dimensional (overall) transverse space, and then $`(ii)`$ compactify the extra (overall) longitudinal directions (if any) and all the relative transverse directions (if any), and perform the Scherk-Schwarz dimensional reduction along the extra (overall) transverse directions (if any). The possible cases are as follows:

* $`\mathrm{\Delta }=4`$ case: M5-brane; NS5-brane; D$`p`$-branes with $`3p8`$.
* $`\mathrm{\Delta }=2`$ case: $`(3|\mathrm{M5},\mathrm{M5})`$; $`(3|\mathrm{NS5},\mathrm{NS5})`$; $`(p1|\mathrm{NS5},\mathrm{D}p)`$ with $`p=4,5,6`$; $`(p2|\mathrm{D}p,\mathrm{D}p)`$ with $`p=5,6`$; $`(p1|\mathrm{D}p,\mathrm{D}(p+2))`$ with $`p=4,5`$; $`(p|\mathrm{D}p,\mathrm{D}(p+4))`$ with $`p=3,4`$.
* $`\mathrm{\Delta }=4/3`$ case: M5$``$M5$``$M5 with 3-dimensional overall longitudinal space; D4$``$NS5$``$NS5; D5$``$D5$``$NS5; NS5$``$NS5$``$D5.

It is interesting to note that the spacetime metric for the 5-dimensional domain wall solution studied in Ref.  is the $`\mathrm{\Delta }=4/3`$ case of the dilatonic domain wall metric. So, the property of the graviton KK modes in such domain wall bulk background is the $`\mathrm{\Delta }=4/3`$ case of what we have studied.

Second, we discuss the string theory embeddings of the 5-dimensional dilatonic domain walls with $`\mathrm{\Delta }4/N`$. To obtain such domain walls, we $`(i)`$ start with (intersecting) branes in ten or eleven dimensions (with equal constituent brane charges) with at least 3-dimensional (overall) longitudinal space and at least 3-dimensional (overall) transverse space, $`(ii)`$ compactify the relative transverse directions (if any), the extra longitudinal directions and possibly some of overall transverse directions to obtain the dilatonic 3-brane in $`D^{}7`$ dimensions, and then $`(iii)`$ compactify on $`S^{D^{}5}`$. The values of the dilaton coupling parameter $`a`$ and $`\mathrm{\Delta }`$ of the resulting 5-dimensional dilatonic domain wall are the $`p=3`$ case of Eq. (14):

$`|a|`$ $`=`$ $`{\displaystyle \frac{2}{3}}\sqrt{{\displaystyle \frac{2(D^{}2)4N(D^{}6)}{2(D^{}5)N(D^{}6)}}},`$ (15)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{4(D^{}6)}{2(D^{}5)N(D^{}6)}}.`$ (16)
A particularly interesting case is the non-dilatonic domain wall ($`a=0`$ case) of the RS model . It is straightforward to check that $`a`$ in Eq. (16) can be zero only for the $`D^{}=10`$ and $`N=1`$ case, i.e., the $`S^5`$-reduction of D3-brane. This possibility of embedding the RS model was later studied in Ref. . The remaining cases of string theory embeddings of dilatonic domain walls with $`\mathrm{\Delta }4/N`$ are as follows:

* $`N=1`$ case:
(1) starting with M5-brane, one compactifies two of the longitudinal directions and then compactifies the transverse space on $`T^n\times S^{4n}`$ ($`n=0,1,2`$) to obtain 5-dimensional domain wall with $`\mathrm{\Delta }=\frac{4(3n)}{5n}=\frac{12}{5},2,\frac{4}{3}`$;
(2) starting with NS5-brane, one compactifies two of the longitudinal directions and then compactifies the transverse space on $`T^n\times S^{3n}`$ ($`n=0,1`$) to obtain 5-dimensional domain wall with $`\mathrm{\Delta }=\frac{4(2n)}{4n}=2,\frac{4}{3}`$;
(3) starting with D$`p`$-brane with $`3p6`$, one compactifies $`p3`$ of the longitudinal directions and then compactifies the transverse space on $`T^n\times S^{8pn}`$ ($`n6p`$) to obtain 5-dimensional domain wall with $`\mathrm{\Delta }=\frac{4(7pn)}{9pn}=\frac{8}{3},\frac{12}{5},2,\frac{4}{3}`$.
* $`N=2`$ case:
(1) starting with $`(3|\mathrm{M5},\mathrm{M5})`$, one compactifies the relative transverse directions and then compactifies the overall transverse space on $`S^2`$ to obtain 5-dimensional domain wall with $`\mathrm{\Delta }=2`$;
(2) starting with $`(p1|\mathrm{NS5},\mathrm{D}p)`$ with $`p=4,5`$, one compactifies the relative transverse directions and $`p4`$ of the overall transverse directions and then compactifies the overall transverse space on $`S^2`$ obtain 5-dimensional domain wall with $`\mathrm{\Delta }=2`$.


text:
# Sources of CMB Spots in Closed Hyperbolic Universes

## I  <br>INTRODUCTION

Some years ago Cayón and Smoot ; hereafter CS, identified a number of ‘cold’ and ‘hot’ spots in the COBE maps of the cosmic microwave background (CMB), as patches of physical density fluctuations (rather than noise) on the surface of last scattering (SLS). A cold (hot) spot, interpreted as gravitational Sachs-Wolfe effect plus fluctuation of radiation temperature on the SLS, corresponds to an increase (decrease) of matter density . These results were used by the author ; henceforth Paper I, in connection with the possibility of the universe’s spatial section being a closed hyperbolic 3-manifold (CHM).

Tables 1a-b list the galactic coordinates in degrees $`(180<l180,90b90)`$ for the positions of the six overerdense $`(C1C6)`$ and the eight underdense $`(H1H8)`$ CS spots. Given the comoving nature of cosmic geometry, it was argued that those spots, for example the overdense ones, might have evolved into galaxy superclusters which are - or may become - observable in our epoch. The underdense spots would have evolved into the relative voids in the observed structure of the large-scale matter distribution.

| Table 1a. Cold spots | | | |
| --- | --- | --- | --- |
| Spot |  | $`l`$ | $`b`$ |
| $`C1`$ |  | $`99.0`$ | $`57.0`$ |
| $`C2`$ |  | $`85.0`$ | $`40.0`$ |
| $`C3`$ |  | $`21.0`$ | $`45.0`$ |
| $`C4`$ |  | $`73.0`$ | $`29.0`$ |
| $`C5`$ |  | $`81.5`$ | $`59.0`$ |
| $`C6`$ |  | $`86.0`$ | $`33.0`$ |

| Table 1b. Hot spots | | | |
| --- | --- | --- | --- |
| Spot |  | $`l`$ | $`b`$ |
| $`H1`$ |  | $`24.0`$ | $`51.0`$ |
| $`H2`$ |  | $`46.0`$ | $`32.0`$ |
| $`H3`$ |  | $`96.0`$ | $`28.0`$ |
| $`H4`$ |  | $`120.0`$ | $`36.6`$ |
| $`H5`$ |  | $`172.0`$ | $`29.0`$ |
| $`H6`$ |  | $`81.0`$ | $`33.0`$ |
| $`H7`$ |  | $`141.0`$ | $`73.5`$ |
| $`H8`$ |  | $`82.0`$ | $`58.0`$ |

Since the CMB spots are the result of using the COBE satellite’s wide-angle detectors, they represent a smearing of density inhomogeneities on the scale of about ten arc degrees. But we may expect that the small angle detectors in the planned projects MAP and Planck will produce spots that cluster in a way consistent with COBE’s maps. In the bottom-up model of structure formation  these fine-grained overdense spots would evolve into galaxies, which would later produce clusters and superclusters. Therefore the original idea that COBE spots become superclusters (or large-scale voids) seems basically sound.

The purpose of Paper I was to fit a number of CHM’s to the CS spots. Since these are interpreted as density inhomogeneities in the fundamental polyhedron (FP) for the manifold, we hoped the positions of the latter, when compared with those of observed structures and voids, might favor some of those CHM’s as the real cosmic space. This task is of course far from simple: there are infinitely many CHM’s, Earth might be located anywhere inside them, with any orientation with respect to the real sky, not to mention the uncertainties about $`\mathrm{\Omega }_0`$, $`H_0`$, and now also $`\mathrm{\Lambda }`$! But we may expect this situation to improve: the parameters’ ranges will be narrowed, and the several proposals for topology determination (cf. Weeks’s  and Roukema’s  talks in this Session) will eventually indicate what manifolds are the most likely candidates for representing the cosmos.

## II THE SEARCH FOR SEEDS IN THE FUNDAMENTAL POLYHEDRON

In Paper I a lopsided method was adopted to choose a possible source for a given CS spot. The sources shown in Fig. 3 and Tables 2 and 3 of Paper I are usually concentrated in narrow bands of galactic latitude $`b`$. Another problem with Paper I was that we looked for sources inside the maximum injectivity FP, with basepoint at the center (which is the standard in SnapPea ) and the observer displaced from the center; this produces an asymmetry in the distances of the found candidate sources.

Here I make two improvements on this research, both with the help of our geometer guru, Jeff Weeks. One of them is that now $`FP`$ is chosen with basepoint on the observer’s position, which makes the source distribution centered on the the observer. The other one is in the search procedure, which will now be discussed.

Fig. 1 is a crude sketch of $`FP`$ for a CHM isometric to a quotient space $`H^3/\mathrm{\Gamma }`$ as usual, with the basepoint and observer at $`FP`$’s center $`O`$, and the spherical $`SLS`$.

With a point $`P^{}`$ (the center of a spot) $``$ $`SLS`$ we want to find the pre-image $`P`$ of $`P^{}`$ in $`FP`$, that is, $`P^{}=\gamma P`$ for some $`\gamma \mathrm{\Gamma }`$. Let $`\{g_k,k=1n\}`$ be the set of face-pairing generators corresponding to the $`n`$-face $`FP.`$ Points $`P`$ and $`P^{}`$, and cell $`FP^{}=\gamma (FP)`$ are also shown in Fig. 1.

The computer search for $`P`$ and $`\gamma `$ is shown in Fig. 2 as a flow chart, where $`d(A,B)`$ means comoving distance between points $`A`$ and $`B`$; $`g`$ is represented by G, $`g_k`$ by Gk, and inv(G) is the inverse of $`g`$. Referring to this figure, note that (i) if the result is $`\gamma =g=1`$, then $`P=P^{}FP`$, which means that $`FP`$ is large enough to intersect $`SLS`$; and (ii) the calculation of $`\gamma `$ is not necessary, but is useful for checking: it must be $`P^{}=\gamma P`$.

The search for sources was done in Paper I for the ten smallest known CHM’s, numbered HW1-10; they are described in . For each of these the observer’s position and space orientation were chosen randomly, but, as said above, not so randomly distributed sources were found. Here we get better results, as shown for manifold $`m007(+3,1)`$ or HW3 in Table 2. The observer is at position $`(0.3,0.0,0.0)`$ in Klein coordinates relative to the standard $`FP`$ in SnapPea, and axes rotated by $`(4.5669,0.1078,5.3451)`$ in Euler angles in radians, with respect to the axes obtained for the displaced basepoint. (It is interesting to observe that the standard $`FP`$ has 20 faces, against 24 faces in our example.) Table 2 shows $`(Z,l,b)`$ and the path $`g=\gamma ^1`$ from the spot to its source, as a word in the generators ($`g_1g_{24}`$ are here labeled $`G00G23`$, as in the computer printout). We adopt the values $`\mathrm{\Omega }_0=0.3`$ and $`Z_{SLS}=1300.`$

The idea to pursue, which is illustrated in Paper I, is to vary these positions and orientations randomly, until we get sources that match data in catalogs of galaxy superclusters and voids - which, for this to become possible, should reach much deeper space than the present limit of about $`Z=0.12`$ in Einasto et al.’s  list of superclusters.

ACKNOWLEDGEMENTS

I am grateful to Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP-Brazil) for a grant to allow my participation in this Meeting. Also to Jeff Weeks for his continued assistance on geometry matters, and to Laerte Sodré Jr. for a conversation on superclusters.


text:
# Language identification of controlled systems: Modelling, control and anomaly detection

## 1 Introduction

Formal language theory has been used in the past to study autonomous systems, the grammatical complexity of the sequences generated according to some coding being used to characterize the complexity of the dynamical system. References  to  are a representative sample of some of this work. However, most systems used in technological applications are not autonomous, but controlled dynamical systems. Therefore one has to distinguish, from the start, between the information that is generated by the dynamics of the system and the one that depends on the free will of the controlling operator. This has some implications on the nature of the formal entities that are used to code the system variables and on the nature of the languages that are generated. Another important requirement for practical applications is the possibility to extract the grammars by on line learning from the data generated by the system, independently of whether or not an analytic model is known.

In Section 2, after a short introduction to some formal language concepts, we show how the modelling framework for controlled dynamical systems leads, almost uniquely, to a context-dependent grammatical formulation. A learning algorithm is then proposed for on-line generation of the productions of the grammar. The following sections illustrate how this formulation may be used for modelling, control and anomaly detection.

Most of our applications concern electromechanical drives. Using grammatical interpolation, as described in Section 3, one may obtain quite accurate results, even with small grammars. However the main capabilities of a language-based formulation lies in the pattern detection domain. Therefore early detection of anomalies or faulty behaviour is probably the most promising domain of application for this technique.

## 2 Grammars for controlled dynamical systems

### 2.1 Formal language concepts

Many fundamental contributions to the formal language field originate from the work of Chomsky  \- , whose theory of formal grammars had a major influence in the development of the subject. Grammatical inference, that is, the development of algorithms which extract grammatical information from examples, is a concept that goes back to Gold’s work . Since then a great deal of work has been done, which can be found in several surveys  \-

To apply grammatical inference procedures, a dynamical system must be considered as an entity (linguistic source) capable of generating a specific language. The grammar $`G`$ of the language is the set of rules that specifies all the words in the language and their relationships. Once the grammar is found, the grammar itself is a model for the source. To define a grammar $`G`$ one specifies a terminal alphabet, a non-terminal alphabet, a start symbol and a set of productions

$$G=\{\mathrm{\Sigma }_T,\mathrm{\Sigma }_N,S,P\}$$
(1)
Terminal alphabet ($`\mathrm{\Sigma }_T`$): set of symbols that make up the words, a word being a string of symbols.

Non-terminal alphabet ($`\mathrm{\Sigma }_N`$): set of auxiliary symbols that are used to generate the words by the production rules.

Start symbol ($`S`$): a special non-terminal symbol used to start the generation of words.

Productions ($`P`$): set of substitution rules (denoted $`ab`$) used to generate the allowed words in the language.

Example: the grammar

$$\left\{\mathrm{\Sigma }_T=(a,b),\mathrm{\Sigma }_N=\left(A\right),S,P=\left(SbA,AaA,Aa\right)\right\}$$
(2)
generates the language with words consisting of a $`b`$ symbol followed by any number of $`a`$ symbols.

The set of productions encodes the dynamics of the system that generates the language. Any word that can be derived from the start symbol by a sequence of productions of the grammar is said to be in the language generated by the dynamical system.

Grammatical inference is an algorithm by which a grammar is inferred from a set of sample words produced by the dynamical system considered as the linguistic source. Notice that there is not a unique relation between a language and the grammar, because distinct grammars may generate the same language.

Grammatical inference, in general, is the identification of a grammar from a set of positive and negative examples. One may also consider a scheme where a ”teacher” answers questions concerning the language to be inferred. However, in dynamical system identification by grammatical inference, we are in the restricted setting of language identification from positive examples alone. Therefore a finite sample cannot characterize an infinite language. It may at most specify the class of languages, which possesses that finite sample as an allowed one.

A practical problem, when learning from a finite set of examples, is the fact that some relevant productions may never be inferred. If a metric is defined in the space of words, this problem may be solved to some extent by grammatical interpolation techniques. Another practical problem arises when the system under study is perturbed by noise. Then, productions may be inferred which are not characteristic of the system, but are a result of the noise perturbation. The solution used in our applications is to keep as valid only those productions that appear many times in the learning process. The threshold for acceptance of a production depends of course on the size of the training set and the estimated noise intensity.

We work here in the phrase-structure grammar setting as defined by Chomsky. Notice however that more general systems like programmed grammars, cooperating distributed grammars or contextual grammars  may also find useful applications in the characterization of dynamical systems. In particular, programmed languages and the interpretation of controls as labels, could be a viable alternative to our identification of the control variables with non-terminal symbols.

### 2.2 Grammatical inference in controlled dynamical systems

A model for a controlled dynamical system has the general form

$$\begin{array}{ccc}\stackrel{}{x}\hfill & =\hfill & f(x,U)\hfill \\ y\hfill & =\hfill & h\left(x\right)\hfill \end{array}$$
(3)
or, considering the time divided into discrete steps

$$\begin{array}{ccc}\hfill x_{k+1}& =\hfill & f(x_k,U_k)\hfill \\ \hfill y_k& =\hfill & h\left(x_k\right)\hfill \end{array}$$
(4)
where $`x`$ is the state variable, $`y`$ the output (or observed) variable and $`U`$ the input (or control) variable. Eqs.(3 \- 4) also establish a functional relationship between the output variables at different times

$$\begin{array}{ccc}\hfill y_{k+1}& =\hfill & g(y_k,U_k)\hfill \end{array}$$
(5)
However, in most systems used in technology, not all state variables are observable. Therefore Eq.(5) does not provide a complete specification of the system. In general, specification of the dynamics in terms of the output variables requires a set of functional relationships involving many time steps in the past, namely

$$\begin{array}{ccc}y_{k+1}\hfill & =\hfill & g_0\left(U_k\right)\hfill \\ y_{k+1}\hfill & =\hfill & g_1(y_k,U_k)\hfill \\ y_{k+1}\hfill & =\hfill & g_2(y_k,y_{k1},U_k)\hfill \\ y_{k+1}\hfill & =\hfill & g_3(y_k,y_{k1},y_{k2},U_k)\hfill \\ & \mathrm{}\hfill & \end{array}$$
(6)
It is this structure, which is required by dynamical considerations on actual controlled systems, that leads in a natural way to our proposal of $`p`$type productions.

To develop a grammatical description and a grammatical inference algorithm for controlled dynamical systems three steps are required. First, the quantification of the variables, then the specification of the nature of the productions and finally a learning algorithm to extract the productions from the experimental data.

* Quantification

Quantification refers to the creation of alphabets for the output variable $`y`$ and the control variable $`U`$. In our approach we associate the terminal alphabet $`\mathrm{\Sigma }_T`$ to the output variable $`y`$ and the non-terminal alphabet $`\mathrm{\Sigma }_N`$ to the control variable $`U`$. Let $`n`$ be the number of terminal symbols and $`m`$ the number of non-terminal symbols. A quantification of the variables is made, in a discrete way, dividing the variables range in equal intervals and associating each interval to a symbol in the alphabet

$$\begin{array}{ccccc}y\hfill & \hfill & y_i\mathrm{\Sigma }_T\hfill & & i=1\mathrm{}n\hfill \\ U\hfill & \hfill & U_j\mathrm{\Sigma }_N\hfill & & j=1\mathrm{}m\hfill \end{array}$$
(7)
* Productions

$`p`$type productions are defined to be substitution rules of the form

$$y_1\mathrm{}y_pU_ky_1\mathrm{}y_py_{p+1}\delta $$
(8)
where $`y_1\mathrm{}y_p`$ is a sequence of terminal symbols, $`U_k`$ a non terminal symbol and $`\delta `$ a special non terminal. $`\delta `$ is used to allow the conclusion, or not, of a generated word, by the use of following special set of productions

| $`\delta `$ | $``$ | $`U_j`$ | , | $`j=1\mathrm{}m`$ |
| --- | --- | --- | --- | --- |
| $`\delta `$ | $``$ | $`\lambda `$ |  |  |

(9)
where $`\lambda `$ denotes the empty symbol.

A $`p`$type production codes the evolution of the output variable, depending on its $`p`$ past values and on the value of the control variable $`U`$. There is therefore a functional relationships between the dynamics of the system and the $`p`$type productions.

* Learning algorithm

To obtain a sample of the language, a sequence of control signals is applied to the system is such a way that the output variable $`y`$ takes values in a sufficiently wide region. The signal evolution is then quantified as described above and the following learning algorithm is used:

1. A $`0`$type production is assumed for every newly occurring control symbol.

2. A new $`(n+1)`$type production is generated each time the data conflicts with the previously established $`n`$type productions. The conflicting $`n`$type productions are also promoted to $`\left(n+1\right)`$type productions or are deleted if there is not sufficient information on the past to do so.

To be able to revise the conflicting $`n`$type productions, promoting them to $`\left(n+1\right)`$type productions, we need to keep in a memory window the record of a certain numbers of past steps in the system evolution. The length of the memory window need not be larger than $`p_{\mathrm{max}}`$ ($`p_{\mathrm{max}}`$ being the highest order allowed for the productions). In general there is not great harm in having zero-length for the memory window and just deleting all the conflicting lower order productions. This is because all the most relevant productions will reappear some time again in the future.

As an example, consider the following symbol sequence obtained from a data sample

$$\begin{array}{ccccccccc}U\hfill & \text{variable}\hfill & :\hfill & A\hfill & B\hfill & A\hfill & A\hfill & B\hfill & A\hfill \\ y\hfill & \text{variable}\hfill & :\hfill & e\hfill & d\hfill & c\hfill & b\hfill & d\hfill & e\hfill \end{array}$$
(10)
At time zero the algorithm analyzes the leading symbols of both the control and the output variable. Since no other information is yet available, a $`0`$type production $`Ae\delta `$ is assumed. After analyzing the second symbol, the algorithm establishes another $`0`$type production $`Bd\delta `$. The third symbol would yield a $`0`$type production $`Ac\delta `$. However, this production contradicts the previously established $`Ae\delta `$, because the grammar is not stochastic. Therefore, in this case, a $`1`$type production $`dAdc\delta `$ is obtained. The $`0`$type production $`Ae\delta `$ is deleted because no information is available on the past of the first symbol ($`e`$). At the next step we obtain the production. $`cAcb\delta `$. Another conflict arises when we reach the last symbol in (10). Then a $`2`$type production $`bdAbde\delta `$ is obtained and $`dAdc\delta `$ is revised to $`edAedc\delta `$. Finally one is left with

$$\begin{array}{ccc}\hfill B& \hfill & d\delta \hfill \\ \hfill cA& \hfill & cb\delta \hfill \\ \hfill bdA& \hfill & bde\delta \hfill \\ \hfill edA& \hfill & edc\delta \hfill \end{array}$$
(11)
The learning algorithm is described by the flow chart in Fig.1.

Unlike, for example, the rules for fuzzy logic inference, the structure of the formal language productions, obtained by our algorithm, is not established in advance. A mixture of different types of productions may be obtained in the grammar. This feature provides flexible adaptation in dynamical systems displaying distinct behavior in different regions of the working space.

In a controlled system the two types of variables, output variables and control variables, have a distinct nature. This is represented in our formalism by the assignment of terminal symbols to the dynamical variables and non-terminal symbols to the control. In turn, the productions represent the action of the control in the context of the past and present dynamics of the system. The languages generated in this way fall in the category of context-sensitive languages.

In practical applications the grammatical inference algorithm is completed with grammatical interpolation and a noise rejection mechanism. Grammatical interpolation is discussed in Section 3 and noise rejection is simply implemented by keeping only those productions that appear a sufficient number of times in the data. Also, for practical purposes, an upper bound is put on the order of the productions that the algorithm generates.

As a simple example, of how this methodology is used to model the dynamics, we consider the following piecewise linear system

$$y_{k+1}=y_ky_{k1}+2U_k\text{ with }\{\begin{array}{c}y_k\{1,2,3,4\}\\ U_k\{1,0,1\}\end{array}$$
(12)
The working domain of this dynamical systems is spanned by three distinct planes, depending on the values of the $`U_k`$ variable (Fig.2).

The assumed codification is

| $`y_k\mathrm{\Sigma }_T`$ | $`U_k\mathrm{\Sigma }_N`$ |
| --- | --- |
| 1 $``$ a | $``$1 $``$ A |
| 2 $``$ b | 0 $``$ B |
| 3 $``$ c | 1 $``$ C |
| 4 $``$ d |  |

The language generated by this dynamical system is represented by a grammar $`G_1`$, whose rewriting system has one $`0`$type production, two $`1`$type productions and twenty four $`2`$type productions.

This system was used to test our on-line learning algorithm with a random input control signal $`U_k`$. Fig.3 shows the evolution of the productions, extracted from randomly generated data, for the first 30 data points. One sees, for example, at the data point 19, the cancellation of a $`0`$type production forced by a new $`1`$type production. Also at data point 62 two contradictory $`2`$type productions are discarded.

## 3 Modelling an electrical drive

The linguistic coding described before has been used, in our laboratory, to model an experimental drive system, consisting of an electronically fed induction machine.

The experimental system is depicted in Fig.4. It is composed of an induction motor driven by a power inverter. Drive systems of this type are modelled using the electromechanical power conversion theory . A 12-equations model represents the electrical drive system. The large number of variables hinders any attempt to perform on-line learning without huge computational costs. Therefore, one needs to simplify the electrical drive model.

Considering an input variable internal control loop , stator currents may be assumed to be controlled. The model is then reduced from a 12th to a 3rd order system  which, in the rotor flux reference frame, is

$$\begin{array}{ccc}\underset{dr}{\overset{}{\psi }}\hfill & =\hfill & \frac{1}{\tau _r}\psi _{dr}+\left(\omega _R\omega \right)\psi _{qr}+\frac{M}{\tau _r}i_{ds}\hfill \\ \underset{qr}{\overset{}{\psi }}\hfill & =\hfill & \left(\omega _R\omega \right)\psi _{dr}\frac{1}{\tau _r}\psi _{qr}+\frac{M}{\tau _r}i_{qs}\hfill \\ \stackrel{}{\omega }\hfill & =\hfill & \frac{M}{JL_r}\left(i_{qs}\psi _{dr}i_{ds}\psi _{qr}\right)\frac{B}{J}\omega T_{ext}\hfill \end{array}$$
(13)
The rotor fluxes ($`\psi _{dr}`$,$`\psi _{qr}`$) and the rotor speed ($`\omega `$) are the state variables, and the stator current components ($`i_{ds}`$,$`i_{qs}`$) the control variables. $`\tau _r`$ is the rotor time constant, $`M`$ the mutual (stator-rotor) induction coefficient, $`J`$ the inertia coefficient, $`L_r`$ the rotor self-inductance coefficient, $`B`$ the friction coefficient and $`T_{ext}`$ the load torque. With the reference frame rotor flux given by

$$\omega _R=\omega +\frac{M}{\tau _r}\frac{i_{qs}}{\psi _r}$$
(14)
one obtains from (13)

$$\begin{array}{ccc}\underset{qr}{\overset{}{\psi }}\hfill & =\hfill & \frac{1}{\tau _r}\psi _r+\frac{M}{\tau _r}i_{ds}\hfill \\ \stackrel{}{\omega }\hfill & =\hfill & \frac{M}{JL_r}i_{qs}\psi _r\frac{B}{J}\omega T_{ext}\hfill \end{array}$$

To obtain relevant training sets, two types of information, qualitative and quantitative, should be acquired. Qualitative information is obtained from the mathematical models. This leads to the choice of the relevant variables that describe the system behavior. This choice is very important to insure the existence of a functional relationship representative of the drive behavior. Otherwise both learning and recognition would be impossible.

To obtain quantitative information experimental data is acquired. To obtain the data an excitation signal must be chosen. A possible approach would be to use a pseudo-random binary signal. However, this signal is not the best choice for drive systems because it is filtered by mechanical time constants. It is better to use sinusoidal signals of different amplitudes and frequencies. In this way, frequencies and amplitudes may be adjusted within the limits of the drive response, thus avoiding the filtering problem and collecting data that adequately spans the operating domain.

A grammar of the drive language is inferred from the (control) input - output experimental information. The $`i_{ds}`$ current is kept constant, a non terminal alphabet is established from the quantification of the $`i_{qs}`$ current and a terminal alphabet established from the output speed signal ($`\omega `$). A representative training set is obtained with a sinusoidal reference signal containing a combination of different amplitudes and frequencies. Fig.5 shows the evolution, in the input/output space, of the training and test data sets, the test data being displayed in bold.

The acquisition of the training set is an essential step in obtaining a good knowledge of the system to be modelled. The training set must cover a representative part of the entire working domain.

### 3.1 Recognition results

To test the validity of the recognition process, two alphabets were considered. A quantification is established for both input and output variables, yielding a 60 symbol alphabet.

The inferred grammar contains 183 productions distributed in the following way: 5 $`0`$type productions, 156 $`1`$type productions and 22 $`2`$type productions. To avoid the influence of noise on the inferred productions, a production is considered valid only if it appears more than a certain number of times, this parameter depending on the size the training set. This procedure works as a filter on noise perturbations.

Applying this grammar to the test set, we obtain the recognition results shown in Fig.6. The actual electromechanical drive speed is shown as a dotted line and the grammar response as a continuous one.

Apart from the quantification error, the recognition results may be considered satisfactory, because the evolution of the drive speed and the grammar results are similar. However, we must point out that there are some situations, indicated by the white arrows in Fig.6, where the grammar does not provide any answer. This happens when no production was inferred representing that particular input/output relationship. As explained later, a generalization method is used to establish non-learned productions.

For the second alphabet, a quantification interval ten times higher is assumed for both input and output variables, yielding now a 6-symbol alphabet, instead of the previous 60 symbols. The inferred grammar contains 36 productions, distributed in the following way: 0 $`0`$type productions, 4 $`1`$type productions and 32 $`2`$type productions. Applying this grammar to the test set the recognition results are those shown in Fig.7. As before, the electromechanical drive speed is shown as a dotted line and the grammar response as a continuous one. The white arrows denote the absence of applicable productions.

As expected, the quantitative modelling performance deteriorates, there being a larger difference between the measured drive speed and the grammatical reconstruction. However, if the important feature is the qualitative recognition of the drive speed, the results may still be considered as satisfactory. For a qualitative recognition process a limited alphabet is not a serious drawback. However, as in the larger grammar, the problem of non-existent productions is present here.

The discrete quantification performed within the grammatical inference algorithm implies a discrete response when recognition of the drive is performed. In order to improve the quantitative recognition capabilities of the algorithm two solutions may be used. The first is to increase the number of symbols in the alphabet to make a more accurate quantification. Alternatively one may use numerical interpolation of the productions. The first solution increases the number of productions thus slowing the recognition process. The second solution does not increase the number of productions.

To illustrate the operation of the interpolation process consider again the 60-symbol alphabet and another training set which is half the size of the one used for the results in Fig.6. After inferring the correspondent grammar, the recognition results are shown in Fig.8. It is seen that the recognition process has considerable faults because, for the reduced training set, the number of productions that is obtained is considerably smaller. A smaller training set led to a faulty grammar, which fails to recognize some words. However because the symbols in our alphabets are in correspondence with numerical values of the variables domain, a metric may be defined in the space of words. New productions may then be obtained, from those inferred from the smaller sample, by interpolation.

Grammatical interpolation establishes new productions by a structural matching procedure. The main idea of structural matching is based on a measure of similarity between the unknown input pattern and the available data structures. The measure of similarity is the distance between the non-existent production and the nearest productions. Several methods for structural word matching have been reported in the literature. The basic algorithm states that the distance between two words is related to the sequence of edit operations (substitution, insertion, and deletion) required to transform one word into another. For any sequence of edit operations a cost function $`c\left(s\right)`$ is considered

$$c\left(s\right)=\underset{i=1}{\overset{n}{}}c\left(e_i\right)$$

$`c\left(s\right)`$ denotes the cost of a particular sequence $`s`$, and $`c\left(e_i\right)`$ the cost of a particular edit operation. The distance between two words $`x`$ and $`y`$ is defined as the minimum cost of transforming a word into another

$$d(x,y)=\mathrm{min}\left\{c\left(s\right)|s\text{ being a sequence of edit operations transforming }x\text{ into }y\right\}$$

When a word cannot be recognized because there is no production generating some symbol in that word, a grammatical interpolation formula is used to obtain the last terminal symbol $`y_{p+1}`$ in the right-hand side of the new production

$$y_{p+1}=\text{quantification}\left(\left(\underset{i=1}{\overset{n}{}}d_i\right)\underset{i=1}{\overset{n}{}}\left(\frac{1}{d_i}y_{p+1}|_{\text{prod }i}\right)\right)$$

This interpolation formula establishes a weighed average of the available productions of the required order, on an average distance between similar productions. The distance $`d_i`$ is the distance between the words $`y_1\mathrm{}y_pU_k`$ in the left-hand side of the existent productions and the word of the new production.

Fig.9 shows the recognition results for the same grammar as considered in Fig.8, when the grammatical interpolation procedure is applied. The effect of the non-existent productions is reduced, and the number of symbols that are not recognized much smaller. Grey arrows denote the interpolation of non-existent productions and white arrows productions that could not be obtained by interpolation.

By generalization we mean the ability of an algorithms to perform even under working conditions distinct of those observed during the learning phase. The interpolation process as described above provides generalization abilities to the algorithm. However this type of generalization is a local effect in the working space. For example, consider the training set shown in Fig.10.

In this case the training set contains only positive speed values. Fig.11 shows the modelling results that are obtained using this sample data. Recognition fails for some of the domain areas not covered by the training data set. Namely symbols that code for negative drive speed are not recognized. Generalization has only a local effect. Therefore, experimental training data must always cover a significant part of the working domain.

## 4 Anomaly detection

Any grammar $`G`$ codes for the class of patterns that belong to the language $`L(G)`$ the grammar generates. In this way one may use the grammar to recognize well formed strings and reject anomalous ones. Once one has learned the grammar of the strings generated by the system under normal conditions, if at a later time there is some fault in the dynamic system, strings that will be generated are not compatible with the learned grammar. It is also reasonable to assume that this might occur already at an early stage of the anomaly. The language algorithm might therefore be valuable as a tool for early fault detection.

Detection of grammar anomalies is based on the same distance and similarity measures as considered before, for grammar interpolation. The anomaly detection algorithm (Fig.12) computes the distance between the string generated by the dynamical system and the one generated by the learned grammar. If this distance exceeds some threshold, a fault is reported.

As a first example, let us consider the dynamical system represented by Eq.(12), with the same codification and inferred grammar as before. For the anomaly detection test we consider a random input during 100 time steps. A sustained anomaly was simulated between time steps 30 to 60, and a spurious one between time steps 80 to 85. The anomaly simulation is simply the replacement of the parameter 2 by 3 in Eq.(12). The dynamical system changes and therefore some of the strings that are generated do not match the learned grammar. Whenever the distance between the system-generated strings and those of the grammar exceeds a threshold, an anomaly is reported.

Figs.13-15 display the anomaly detection results. In response to the control variable both the learned grammar and the system produce terminal words. These are shown in Fig.14, the output of the grammar as a continuous line and the output of the dynamical system as a dashed line. They differ whenever there is an anomaly. Fig.15 displays the distance between words produced by the learned grammar and those of the dynamical system, for a word length of 10 symbols. This distance exceeds the threshold only for the sustained anomaly.

As a second example we consider anomaly detection in an induction motor. The anomaly is a rotor broken bar, and the variable coded as a terminal alphabet is the stator phase current. The learned grammar, representing the experimental drive system, is generated with the motor operating without anomaly. Fig16 shows the stator current deviation due to the presence of the anomaly between time 1 and 2 (seconds). Fig.17 shows the distance between the words generated by the learned grammar and those of the induction motor. Soon after the onset of the anomaly, this distance begins to exceed the threshold. A non-zero threshold must be defined because, as is seen in the figure, some deviations are obtained even without anomalies due to noise or even to small mismatches in the learned grammar.

## 5 Control

Once the grammar generated by a dynamical system is learned, this grammar, defining the structural features of the words produced by the linguistic source, is a model for the source. So it is possible to use the productions to predict the evolution of the system or to find the control string that leads from one state to another.

The basic idea of the control algorithm is as follows. A string of terminal symbols is given which represents the sequence of states that one wants the system to follow. Knowing the initial state of the system, for each symbol of the desired sequence of target states one chooses a production leading the system from its present state to the desired one. At each step the production that is chosen is the one that leads the system as close as possible to the desired state. The procedure is repeated until that symbol is reached and then one moves to the next target symbol.

As an example consider the dynamical system defined by Eq.(12). Suppose that the sequence of target states is ’$`cbcbaabcccbacc`$’. Fig.18 shows the evolution of the system obtained by the control methodology described above. The target sequence is shown as a dashed line and the system evolution as a continuous line.

Notice that to evolve from one target symbol to another the system may have to pass through other terminal symbols. The control algorithm is robust in the sense that whenever, at some time $`t`$, the expected evolution differs from the one predicted by the learned grammar, due to some perturbation or noise, it suffices to apply the same algorithm once more starting from the state at time $`t`$.

## 6 Conclusions and Remarks

This paper proposes a characterization of controlled system dynamics by formal language techniques. The applications we have been concerned with are mostly electromechanical drives. However, the technique might also be useful for many other dynamical systems or industrial processes.

An important issue is the choice of an appropriate quantification for the feature space. This influences the choice and size of the alphabets and the dimension and complexity of the grammars. Also, good coverage of the working domain, in the learning stage, is essential to insure good generalization properties.

Unlike other approaches, the nature of the rules that define the dynamical system are not set up in advance, different types of productions being established on-line according to the incoming words from the linguist source. The method seems to provide good recognition results and a reasonable quantitative accuracy. Fuzzy logics approaches are probably more efficient when quantitative accuracy is very important and a good knowledge already exists concerning the rules of the system. However, if the dynamics is unknown or if pattern recognition is the important issue, the formal language approach seems quite promising.

Fault detection in control systems, industrial plants and power networks is a field where this work is being pursued and compared with immunity-based learning systems.


text:
# Computational Geometry Column 39

## Abstract

The resolution of a decades-old open problem is described: polygonal chains cannot lock in the plane.

A polygonal chain is a connected series of line segments. Chains may be open, or closed to form a polygon. A simple chain is one that does not self-intersect: only segments adjacent in the chain intersect, and then only at their shared endpoint. If the segments of a polygonal chain are viewed as rigid bars, and the vertices as universal joints, natural questions are whether every open chain can be straightened—reconfigured to lie on a straight line—and whether every closed chain can be convexified—reconfigured to form a planar convex polygon. In both cases, the chains are to remain simple throughout the motion. If a chain cannot be so reconfigured, it is called locked.

These questions were raised by several researchers independently since the 1970’s,111 G. Bergman, U. Grenander, W. Lenhart, J. Mitchell, S. Schanuel, and S. Whitesides. It seems to have first appeared in print in 1995: \[LW95\] and \[Kir95, p. 270\].  and were the subject of intense investigation by the late 1990’s. It was first established that chains can lock in three dimensions (3D): both locked open chains, and locked closed but unknotted chains are possible \[CJ98, BDD+99\]. In 4D, neither open nor closed chains can lock \[CO99\]. But aside from some special cases (e.g., star-shaped polygons cannot lock \[ELR+98\], polygonal trees can lock \[BDD+98\]), the problems remained unresolved for chains in 2D.

Connelly, Demaine, and Rote have now settled the questions, establishing that neither open nor closed chains can lock in 2D \[CDR00\]. Their result is even more general: no collection of disjoint simple chains are locked (although of course a chain nested inside a polygon is forever confined). They prove that every such collection has an expansive motion: one during which the distance between every pair of vertices increases or stays the same. An expansive motion automatically maintains simplicity, for the distance between some pair of vertices would have to decrease to reach self-intersection. Their proof first establishes that any configuration of chains has an infinitesimal expansive motion, and then combines these motions into a global expansive motion. Although there are different ways to stitch the infinitesimal motions together, one produces especially natural movements. This motion minimizes the squared lengths of the velocity vectors $`v_i`$ applied to each vertex $`i`$ at each time: minimize $`_iv_i^2`$, subject to constraints maintaining segment lengths and forcing expansivity.222 Their proof employs a further nonquadratic “barrier” penalty term for smoothness.  The objective function is strictly convex and the system can be solved by quadratic programming. Figs. 2 and 2 illustrate the solution so obtained for a collection of three closed chains (triangles) entangled with one open chain.

Let us simplify the discussion to a single open chain and offer a crude sketch of their key argument. View nonadjacent vertices of the chain as connected by a strut, which is permitted to increase in length or stay the same, but never decrease. If this bar-and-strut framework has an infinitesimal motion, then the motion is necessarily expansive. An equilibrium stress in an assignment of weights (stresses) to each bar and strut so that (a) the stress on each strut is nonnegative, and (b) the stress vectors (vectors along the bars/struts scaled by the stresses) are in equilibrium at each vertex. Bar stresses may be positive (compression) or negative (tension). If a framework is “stuck,” or, more formally, “first-order rigid” \[CW96, Whi97\], then it has a nonzero equilibrium stress. So if a framework only has the trivial zero equilibrium stress, then it possesses an infinitesimal expansive motion.

In order to apply a theorem that holds for planar frameworks, every intersection point of the framework is used to divide the bars/struts, which produces a planar framework equivalent to the original in terms of equilibrium stresses. Now the century-old Maxwell-Cremona theorem is applied: if the framework has a nonzero equilibrium stress, then it can be “lifted” to a nonflat polyhedral terrain that projects to the framework, with positive-stress edges (always bars) lifting to “mountains,” and negative-stress edges to “valleys.” Finally it is shown that such a lifting is impossible, by concentrating on its maximum. For example, if this maximum is a single point, then it must lie at the junction of at least three mountain-edges, which violates the fact that every vertex of the chain is incident to at most two bars. A case analysis on the topologically possible maxima proves that only the degenerate flat lifting is possible, which implies that the original framework has only the zero equilibrium stress, which implies that it has an infinitesimal expansive motion, which implies that it has a global expansive motion, i.e., it is not locked.

Building on this work, Streinu has found a way to decompose the expansive global motion for chains of $`n`$ vertices into $`O(n^2)`$ sections, each of which is the motion of a one degree-of-freedom mechanism \[Str00\]. This mechanism is constructed by removing one hull edge from a pseudotriangulation of the chains, a partition by diagonal bars into regions each bounded by three reflex chains (where a single segment counts as a reflex chain). The mechanism is opened following its only free trajectory in configuration space until two adjacent edges align, at which point the pseudotrianguation is revised locally, and the expansion process continued.


text:
# References

LPT Orsay/00-58

SINP/TNP/00-14

hep-ph/0007016

Can R-parity violation explain the LSND data as well?
Asmaa Abada 1 and Gautam Bhattacharyya 2 1Laboratoire de Physique Théorique, Université de Paris XI, Bâtiment 210,
91405 Orsay Cedex, France

2Saha Institute of Nuclear Physics, 1/AF Bidhan Nagar, Calcutta 700064, India
Abstract

> The recent Super-Kamiokande data now admit only one type of mass hierarchy in a framework with three active and one sterile neutrinos. We show that neutrino masses and mixings generated by R-parity-violating couplings, with values within their experimental upper limits, are capable of reproducing this hierarchy, explaining all neutrino data particularly after including the LSND results.
>
> PACS number(s): 12.60.Jv, 14.60.Pq, 11.30.Fs

The vacuum oscillation interpretation of solar neutrino data requires $`\mathrm{\Delta }m_{\mathrm{solar}}^210^{10}\mathrm{eV}^2`$, while the matter enhanced MSW solution prefers the range $`\mathrm{\Delta }m_{\mathrm{solar}}^2(10^510^4)\mathrm{eV}^2`$ . The atmospheric neutrino anomaly can be explained by $`\mathrm{\Delta }m_{\mathrm{atm}}^2(25)\times 10^3\mathrm{eV}^2`$ (with $`\mathrm{sin}^22\theta _{\mathrm{atm}}>0.88`$) . In the standard three neutrino framework, which offers two independent mass differences, it is possible to conceive of mass hierarchies which can explain both solar and atmospheric results. If in addition the LSND data, carrying a positive indication of $`\nu _\mu \nu _e`$ oscillation with $`\mathrm{\Delta }m_{\mathrm{LSND}}^2(0.31)\mathrm{eV}^2`$ , are also sought to be simultaneously explained, then one has to expand the standard three neutrino scenario as it falls short of one independent mass difference. In fact, within this three neutrino framework one can fit any two of the three $`\mathrm{\Delta }m^2`$ mentioned above. Under the circumstances, one often leaves the LSND data out of consideration pending further confirmation from the MiniBooNE  at FNAL or MINOS long baseline  experiments, since KARMEN  does neither confirm nor exclude the LSND results. On the other hand, if one includes the LSND results to be explained together with the solar and atmospheric data, the minimal extension of the standard scenario that needs to be done is to add one sterile neutrino ($`\nu _s`$)  to the list of three active states. The resulting four neutrino picture can explain all the data .

What are the possible choices of mass hierarchies among these four neutrinos? The data suggest that the choices are very limited. First, we devide them into two types, I and II. In type I, there are three almost degenerate states explaining the solar and atmospheric results, with a fourth state separated by a large gap. This hypothesis cannot explain the LSND results, since the LSND data require that the separated fourth state will have to be either $`\nu _e`$ or $`\nu _\mu `$, but each of them will have to be closely spaced with a third state to explain the solar and atmospheric results. If one disregards the LSND results till further confirmation, then of course type I scenario is allowed with the isolated state either $`\nu _s`$ or $`\nu _\tau `$. In type II, there are two pairs of approximately degenerate states separated by a large gap. The mixing between the two pairs is very small. This scenario can explain all the data, or in other words, the large gap could be the LSND gap. Two cases may arise in this framework. In type IIa scenario, ($`\nu _\tau `$-$`\nu _e`$) form one pair which explains the solar neutrino data, while ($`\nu _\mu `$-$`\nu _s`$) form the other pair maximally mixed to explain the atmospheric anomaly. In type IIb, ($`\nu _\tau `$-$`\nu _\mu `$) form the pair that explains the atmospheric anomaly, while ($`\nu _s`$-$`\nu _e`$) pair explains the solar data 111$`\nu _e\nu _s`$ solar neutrino oscillation will be tested by the SNO experiment soon .. The most recent Super-Kamiokande (SK) data rule out the $`\nu _\mu \nu _s`$ interpretation of atmospheric neutrino anomaly at 99% CL , thus strongly disfavouring the type IIa scenario. We are then left with type IIb as the only surviving option, shown in Fig. 1222This framework is consistent with the Big-Bang Nucleosynthesis constraints on sterile - active mixings .. It should be noted that the oscillation data cannot discriminate between cases that occur by interchanging either the members within a given pair or the relative spacing of the two pairs. These cases are not separately shown.

Under the circumstances, it is a timely exercise to identify those scenarios which reproduce the type IIb spectrum. Do the R-parity-violating ($`\overline{)}R`$) supersymmetric models  fall in that category? In this note we seek to find an answer to this question. Defined in terms of lepton number ($`L`$), baryon number ($`B`$) and spin ($`S`$) of the particle as $`R=(1)^{3B+L+2S}`$, this discrete parity is +1 for all SM particles and -1 for their superpartners. Since neither $`L`$\- nor $`B`$-conservation is ensured by gauge invariance or any such fundamental principles, R-parity is an ad hoc symmetry put in by hand. Although there is no experimental confirmation yet in favour of non-vanishing $`\overline{)}R`$ interactions, the neutrino oscillation data are somehow suggestive that it would be premature to abandon those couplings, as the origin of neutrino masses and mixings could be traced to some of these non-vanishing $`L`$-violating couplings. In order not to allow rapid proton decay we do not switch on $`L`$ and $`B`$ violations simultaneously. The following $`L`$-violating terms in the superpotential are then allowed: (a) $`\lambda _{ijk}L_iL_jE_k^c`$, (b) $`\lambda _{ijk}^{}L_iQ_jD_k^c`$, and (c) $`\mu _iL_iH_u`$. In these expressions, $`L_i`$ and $`Q_i`$ are SU(2)-doublet lepton and quark superfields, $`E_i^c`$ and $`D_i^c`$ are SU(2)-singlet charged lepton and down quark superfields, and $`H_u`$ is the Higgs superfield that gives masses to up-type quarks. Stringent experimental constraints on these a priori independent parameters exist in the literature .

Attempts to fit the observed neutrino data by masses and mixing angles generated by $`\overline{)}R`$ couplings have been done in the past, both in the context of bilinear  as well as trilinear  $`L`$-violating parameters. Most of the analyses have been carried out in the three-neutrino picture. A few remarks on those analyses are in order. Most of them discarded the LSND data as something not so reliable, and confined the discussions to the possibility of fitting only solar and atmospheric data. Adhikari and Omanovic , on the other hand, tried to fit solar, atmospheric (preliminary SK) and LSND all together in a three-neutrino picture and that too considering only trilinear $`\overline{)}R`$ couplings. But they did not consider the SK zenith angle dependence and also assumed an energy independent solar neutrino solution by ignoring the Chlorine data of the Homestake mine experiment. In fact, if one takes all present data into consideration, it is not enough just to add bilinear $`\overline{)}R`$ terms to their analysis – the data compel one to go for a four-neutrino picture. Ref.  deals with the simultaneous presence of bilinear and trilinear parameters. Although the emphasis in ref.  is mostly on finding an explanation for the solar and atmospheric neutrino data in a three-neutrino framework, a qualitative discussion on how to simultaneously explain the LSND data by admitting a fourth sterile state has also been presented.

In the present work, we perform a numerical study to examine whether R-parity violation, with bilinear and trilinear parameters together, can reproduce the type IIb spectrum. Even though the data allow an interchange of the location of the two pairs, we work in a situation where ($`\nu _\mu `$-$`\nu _\tau `$) form the heavier pair333This hierarchy is stable under radiative corrections .. We assume that the masses of all active neutrinos are generated by R-parity violation. Since the mixing between the two pairs is small, we can, for all practical purposes, focus on the heavier pair and work in the $`\nu _\mu `$$`\nu _\tau `$ subspace, always assuming that the parameters $`\lambda _{1jk}^{}`$ and $`\mu _1`$, responsible for the $`\nu _e`$ mass generation, are much smaller. We also assume that the sterile state receives mass of the order of the $`\nu _e`$ mass from some other source. More specifically, we concentrate on those parameters which generate the masses and mixing angles in the above ($`2\times 2`$) subsector, keeping in mind that each of the two absolute masses should be of the order of the LSND gap. In other words, we parametrize the $`\nu _\mu `$$`\nu _\tau `$ mass matrix in terms of the $`\overline{)}R`$ couplings, vary them within physical ranges, and then observe whether there exist solutions that simultaneously satisfy the following experimental constraints:

$`\{\begin{array}{cc}\mathrm{\Delta }m_{\mathrm{LSND}}^2m_4^2m_3^2(0.31)\mathrm{eV}^2,\hfill & \\ \mathrm{\Delta }m_{\mathrm{atm}}^2=m_4^2m_3^2(25)\times 10^3\mathrm{eV}^2,\hfill & \\ \mathrm{sin}^22\theta _{\mathrm{atm}}\mathrm{sin}^22\theta _{34}>0.88.\hfill & \end{array}`$ (4)
Now we turn our attention to the analytic expressions of the neutrino masses induced by bilinear and trilinear couplings. We write the mass matrix as $`_\nu =_\nu ^{\mathrm{tree}}+_\nu ^{\mathrm{loop}}`$. The bilinear couplings contribute to the tree mass . In a basis where there are no sneutrino vacuum expectation values (VEVs) , this can be expressed as444In principle, one can rotate away the bilinear $`\mu _i`$ terms from the superpotential, but still having the $`\mathrm{\Delta }L=2`$ effects via the presence of sneutrino VEVs after the minimization of the scalar potential. One can as well work in a framework where there are no sneutrino VEVs but the $`\mu _i`$ parameters are present in the theory . We work in the latter basis. We emphasize that going from one basis to another, the parametrization would change, but the general conclusion we draw at the end remains unaffected. For a discussion on the basis-independent parametrizations of R-parity violation, see refs. .

$`_{\nu _{ii^{}}}^{\mathrm{tree}}=g_2^2{\displaystyle \frac{(M_1+\mathrm{tan}^2\theta _WM_2)}{4detM}}\mu _i\mu _i^{}v_d^2\alpha \mu _i\mu _i^{}(v_d=H_d^0),`$ (5)
where $`M_{1,2}`$ are the gaugino masses, and $`detM`$ is the determinant of the $`(4\times 4)`$ neutralino mass matrix in the R-parity-conserving case. Considering only the $`\lambda ^{}`$ couplings, the one-loop squark-mediated contribution to the neutrino mass, in the same basis as before, can be written as

$$_{\nu _{ii^{}}}^{\mathrm{loop}}\frac{N_c\lambda _{ijk}^{}\lambda _{i^{}kj}^{}}{16\pi ^2}m_{d_j}m_{d_k}\left[\frac{f(m_{d_j}^2/m_{\stackrel{~}{d}_k}^2)}{m_{\stackrel{~}{d}_k}}+\frac{f(m_{d_k}^2/m_{\stackrel{~}{d}_j}^2)}{m_{\stackrel{~}{d}_j}}\right],$$
(6)
where $`f(x)=(x\mathrm{ln}xx+1)/(x1)^2`$. Here, $`m_{d_i}`$ is the down quark mass of the $`i`$th generation, $`m_{\stackrel{~}{d}_i}`$ is an average of $`\stackrel{~}{d}_{Li}`$ and $`\stackrel{~}{d}_{Ri}`$ squark masses, and $`N_c=3`$ is the colour factor. While writing Eq. (6), we assumed that the left-right squark mixing terms are family-diagonal and are proportional to the corresponding quark masses, i.e., $`m_{\mathrm{LR}}^2(i)=m_{d_i}m_{\stackrel{~}{d}_i}`$. The expression of the $`\lambda `$-induced slepton-mediated contributions to the neutrino mass is similar to Eq. (6), and we do not display it here555We neglect the one-loop diagrams induced by the product of bilinear and trilinear couplings ..

In the basis $`(\nu _\mu ,\nu _\tau )`$, the mass matrix can be parametrized as

$`_\nu =\left(\begin{array}{cc}K\lambda _2^2+\alpha \mu _2^2& K\lambda _2\lambda _3+\alpha \mu _2\mu _3\\ K\lambda _2\lambda _3+\alpha \mu _2\mu _3& K\lambda _3^2+\alpha \mu _3^2\end{array}\right),`$ (9)
where $`\alpha `$ is given by Eq. (5), $`\lambda _{2,3}`$ are two generic trilinear couplings, and $`K`$ captures the loop factors that enter into Eq. (6).

We observe that the determinant of the mass matrix with only bilinear or with only trilinear couplings is identically zero666If we have more than one combination of trilinear couplings, then although the determinant of the mass matrix with trilinear couplings alone will not be zero, still much more fine-tuning may be necessary in order to arrive at any solution with only trilinear couplings.. But this yields a big hierarchy between $`m_3`$ and $`m_4`$, contrary to the requirement that these two states should be approximately degenerate. Only after taking the bilinear and trilinear parameters together we obtain solutions that pass the test in Eq. (4). Actually we need at least four independent input parameters in the mass matrix in order to fit the data. In this case, they are the two trilinear couplings ($`\lambda _2`$, $`\lambda _3`$) and the two bilinear mass parameters ($`\mu _2`$, $`\mu _3`$). Notice that these four parameters as well as $`\alpha `$ can take either sign, while $`K`$ is always positive. We point out that the tree and loop contributions may have the same or opposite signs, and this relative sign plays a crucial role in deciding which combinations of parameters are allowed by the data. After diagonalizing the mass matrix in Eq. (9), we demand that the eigenvalues $`m_{3,4}`$ and the mixing angle $`\theta _{34}=\theta _{\mathrm{atm}}`$ satisfy Eq. (4). In Figs. 2a and 2b we have displayed only a part of the solutions by plotting the acceptable mass spectrum as a function of some allowed input parameters, just to demonstrate that the mechanism of R-parity violation works as a viable explanation777In principle, the allowed spectrum should have been displayed as a five-dimensional plot where all the four input parameters are varying. Mostly for the purpose of a simplified presentation, merely to point out that there indeed exist solutions satisfying Eq. (4), we plotted the spectrum in two dimensions in Figs. 2a and 2b. In each plot, the ‘other three’ parameters are also varying.. The above parametrization is rather general as one can apply it for any $`\lambda ^{}`$ (or $`\lambda `$, for that matter) couplings irrespective of the second and third generation indices. Also at this stage one need not specify the squark (or slepton) and gaugino masses as they are absorbed in $`K`$ and $`\alpha `$ respectively. We observe that the ‘filter’ of Eq. (4) prefers a negative $`\alpha `$, which implies that one or more gaugino mass parameters could be negative (see Eq. (5)). We stress though that we do get some solutions with positive values of $`\alpha `$ as well.

Now assuming, as an illustrative example, that $`\lambda _{233}^{}`$ and $`\lambda _{233}^{}`$ are the only dominant trilinear couplings, the factor $`K`$ turns out to be $`KN_cm_b^2/8\pi ^2m_{\stackrel{~}{b}}`$. Taking the squarks and gauginos to be approximately at the $`300`$ GeV scale, one obtains $`\alpha 2.10^4\mathrm{GeV}^1`$, and $`K1.10^3`$ GeV. The minimum and maximum values of the input parameters that pass this test turn out to be

$`\lambda _{233}^{}\lambda _{333}^{}:[1.3\times 10^3,1.3\times 10^3];`$
$`\mu _{2,3}(\mathrm{GeV}):[5.0\times 10^3,1.0\times 10^3].`$ (10)
These values are consistent with the existing constraints on the above parameters .

To conclude: If R-parity violation has to explain all neutrino data, it is essential to have both trilinear and bilinear $`\overline{)}R`$ terms in addition to having a sterile neutrino in the model. Then it is possible to generate maximally mixed $`\nu _\mu `$ and $`\nu _\tau `$ with their absolute masses in the eV range and mass-squared difference in the milli-$`\mathrm{eV}^2`$ range, as required by the data. The quoted ranges of the $`\overline{)}R`$ parameters within which we obtain solutions are based on certain simple-minded but plausible approximations made for the ease of presentation. The bottom line of our analysis is that we provide an affirmative answer to the question we have asked in the title.

Note added: While we were finishing this note, we became aware of a preliminary SK solar neutrino analysis update  disfavouring a pure $`\nu _e\nu _s`$ solar neutrino oscillation at 95% CL. First, this is only a 2-$`\sigma `$ result which is not enough to exclude a model. Second, it has been claimed that a rate $`+`$ spectrum combined analysis in a full four-flavour scenario exhibits an allowed zone . So the fate of the sterile state may not be that dwindling, and we must wait till the SNO experiment tests this option.

Both authors thank the CERN Theory Division, where this work has been done, for its warm hospitality. They also thank H. Dreiner and A. Raychaudhuri for illuminating comments on the manuscript. GB also thanks LPT, Orsay, for its hospitality, and S. Choubey and S. Goswami for clarifying remarks on neutrino data.


text:
# Microscopic dynamics of a phase transition: equilibrium vs out-of-equilibrium regime

## 1 INTRODUCTION

In the last years there has been a lot of interest in the nuclear physics community for multifragmentation reactions and for the study of liquid-gas phase transition in finite systems such as nuclei . It is still a debated question whether nuclear multifragmentation is an equilibrium type of transition, and if it is first or second-order . With the scope of adding new useful arguments to this debate, we present for the first time to the nuclear physics community the Hamiltonian Mean Field (HMF) model, a system of classical spins coupled through long-range attractive forces. HMF can be solved analitically in the canonical ensemble and shows a second-order phase transition in the thermodynamic limit. Moreover HMF has the double advantage of allowing a microcanonical and dynamical approach to explore the dynamics of a phase transition in a finite system. We believe that, in the same spirit of the Ising model, introduced many years ago, but still studied with extreme interest in statistical mechanics, the HMF model, though its dynamics is probably less complex than that one of a real system, can become a tool of primary importance to extract information on the microscopic dynamics of a second-order phase transition.

In this paper we discuss the equilibrium and out-of equilibrium properties of the model, focusing our attention on the importance of the relaxation to equilibrium and on the possible discrepancies between different ensembles. In this respect the model can be extremely useful also for the nuclear multifragmentation phase transition. In particular in section 3, we show how misleading information on the order of a phase transition can be obtained from a specific heat analysis  if the system is not perfectly equilibrated.

## 2 HMF: THE CANONICAL ENSEMBLE

The Hamiltonian Mean Field (HMF) model was introduced by Ruffo  in 1994 and has been since then intensively studied both analitically and numerically . The Hamiltonian of the model describes a system of N fully-coupled particles moving on the unitary circle, where the particles are characterized by the angles $`\theta _i`$ and the conjugate momenta $`p_i`$. It reads:

$$H(\theta ,p)=K+V=\underset{i=1}{\overset{N}{}}\frac{p_{i}^{}{}_{}{}^{2}}{2}+\frac{1}{2N}\underset{i,j=1}{\overset{N}{}}\left[1cos\left(\theta _i\theta _j\right)\right],$$
(1)
being $`K`$ and $`V`$ the kinetic and potential energy. If one considers a spin vector associated to each particle $`𝐦_i[cos\theta _i,sin\theta _i]`$, the Hamiltonian then describes a linear chain of $`N`$ classical fully-coupled spins, similarly to the XY model, and we can define a total magnetization vector $`𝐌[M_x,M_y]=\frac{1}{N}_{i=1}^N𝐦_i`$ .

A very important feature of the HMF model is that it presents an analytical solution in the canonical ensemble. The solution has been worked out by Antoni and Ruffo, by means of the Hubbard-Stratonovich method to perform the integration over the angles. The system behaves as a ferromagnet at low energy and shows a second-order phase transition at the critical energy density $`U_c=0.75`$ ($`U=E/N`$ where $`E`$ is the total energy), corresponding to a critical temperature $`T_c=0.5`$ . If the sign of the interaction is inverted, one gets an antiferromagnetic behavior: this case has been studied in detail in refs . The theoretical caloric curve reads:

$$U=\frac{E}{N}=\frac{(\beta F)}{\beta }=\frac{1}{(2\beta )}+\frac{1}{2}\left(1M^2\right),$$
(2)
where $`\beta =1/kT`$ is the inverse of the temperature (the Boltzmann constant $`k`$ is set to be equal to $`1`$). The canonical caloric curve is shown in fig.1 (a) as a full curve: it increases almost continuously with energy, but has a discontinuity at the critical point. The modulus of the total magnetization, $`M`$, shown as full curve in fig.1(b), is close to one for very low energy and decreases to zero by increasing the energy. It is equal to zero for $`UU_c`$. It has been shown that the phase-transition in the HMF model is a second-order one with mean field critical exponents $`\beta =1/2`$, $`\alpha =0`$ . It is very intriguing the fact that, though HMF was not constructed to describe a nucleus, its caloric curve turns out to be very similar to the one measured by the ALADIN group and more recently by other experimental collaborations for nuclear systems.

## 3 HMF: THE MICROCANONICAL ENSEMBLE

In this section we discuss the microcanonical numerical simulations and we compare with the equilibrium canonical solution. The Hamilton equations of motion for the $`N`$ particles are

$$\dot{\theta _i}=p_i,\dot{p_i}=sin(\theta _i)M_x+cos(\theta _i)M_y,i=1,\mathrm{},N,$$
(3)
and can be solved numerically for different sizes of the system. In this way HMF allows a dynamical microcanonical approach, perfect to explore the dynamics of a phase transition in a finite system. The equations of motion were integrated by means of a 4th order simplectic algorithm  with a time step $`\delta t=0.2`$ and a relative error in the total conserved energy smaller than $`\mathrm{\Delta }E/E=10^5`$. More details can be found in Refs .

### 3.1 Equilibrium Results

In this paper we consider a system with $`N=500`$, a number of particles of the same order of the typical number of nucleons in a nuclear multifragmentation reaction. We start the system out-of-equilibrium, in the so-called “water bag”, i.e. by putting all rotators at $`q_i=0`$ and giving an initial velocity according to a constant probability distribution function of finite width centered around zero . We remind the reader that, though the interaction in HMF is not explicitly constructed to model a nuclear system, such a far-off-equilibrium initial condition is a good way to simulate the strong off-equilibrium conditions present in a hot and compressed nuclear system before multifragmentation.

We follow the dynamics of HMF up to a final time $`t=210^6`$ and we show in fig.1 the final results of the numerical integration. The numerical results were integrated also on several different runs. Such long integration times are necessary because, when out-of-equilibrium initial conditions are used, the system equilibration can be very slow. This is especially evident close to the critical point, where the presence of quasi-stationary (long living) non-equilibrium states (QSS) has been revealed. The temperature, the total magnetization, the kinetic fluctuations and the specific heat per particle are reported as a function of the energy density as open squares and compared to the theoretical predictions (full curves). The temperature is computed from the average kinetic energy per particle $`T=2<K>/N`$, where the symbol $`<>`$ stands for time averages. The kinetic energy fluctuations $`\mathrm{\Sigma }`$, obtained from the scaled variance of the kinetic energy $`\mathrm{\Sigma }=\frac{\sqrt{<K^2><K>^2}}{\sqrt{N}}`$, are compared to the theoretical prediction in the microcanonical ensemble. The latter is given by the formula $`\mathrm{\Sigma }=\frac{T}{\sqrt{2}}\sqrt{1[12M(\frac{dM}{dT})]^1}`$ . Finally, in panel (d), we compute the specific heat per particle at constant volume in the microcanonical ensemble using the formula proposed for the first time by Lebowitz, Percus and Verlet :

$$C_V=\frac{1}{2}\left[12\left(\frac{\mathrm{\Sigma }}{T}\right)^2\right]^1.$$
(4)
In general, the numerical microcanonical caloric curve and the magnetization obtained at very long times agree quite well with the canonical predictions, apart from small finite-size effects located mainly around the critical point and in the supercritical region, the homogeneous phase. In the latter, the magnetization, due to the fact that N=500, is not exactly zero and assumes a finite value, i.e. $`M1/\sqrt{N}`$. This has visible consequences on the temperature, according to formula (2) and thus also on the kinetic energy fluctuations and the specific heat. We note that the microcanonical specific heat at equilibrium is always positive and agrees with the theoretical formula (reported as a solid line). In conclusions the canonical ensemble and the microcanonical simulations at equilibrium are in good agreement. The equivalence of different ensembles in the thermodynamic limit is also supported by recent analytical results in the microcanonical ensemble .

### 3.2 Out-of-Equilibrium results

In the previous subsection we have presented the microcanonical numerical results obtained in the equilibrium regime. In the present one we study the process of relaxation to equilibrium and we show that HMF has a very rich dynamics.

In fig.2 we report the time evolution of the temperature for different energies, and we study how it relaxes to the canonical value, represented by dashed straight lines. Though for U=0.8,0.4 the equilibrium temperature is almost reached already at a time $`t=3000`$, close to the critical point, and in particular for the energy U=0.7 , we have to wait for much longer times (for a time $`t>510^5`$). In fact, before the relaxation to the canonical value, there is a well defined plateau, where the temperature assumes a constant value smaller than the canonical prediction. The same behavior appears in a whole energy region $`0.5<UU_c`$ and, as shown in ref., it corresponds to the presence of Quasi-Stationary States (QSS). The plateaus are longer the greater the system size, and in the $`N\mathrm{}`$ the system does not relax to the standard canonical equilibrium, remaining forever in the QSS . The QSS represent a different dynamical equilibrium for non-extensive systems, and according to an intriguing scenario recently proposed by Tsallis , the reason could be the inversion of the time limit with the size limit, an inversion which is implicit in our numerical simulations.

In fig.3 we report the same quantities shown in fig.1, now calculated for a very short time, i.e. at $`t=3000`$ instead of $`t=210^6`$. For this time, average variables, such as $`T`$ and $`M`$, have almost reached the canonical equilibrium for a large range of initial energies $`U`$, see panel (a) and (b). There are, however, relevant differences with respect to the equilibrium values in the region $`0.5<UU_c`$, in correspondence to the presence of the QSS, where the caloric curve T(U) shows a backbending typical of a first-order phase transition, see in particular the magnification reported in the inset of panel (a). On the other hand, fluctuations are affected by non-equilibrium effects in the whole spectrum of energy. In fig.3(c) the kinetic energy fluctuations are bigger than the equilibrium ones. As a consequence, see eq.4, the specific heat per particle reported in fig.3(d) assumes negative values and the behavior of $`C_V`$ versus U thus obtained is similar to the one found in nuclear multifragmentation data, see fig.4 of ref.. In particular the authors of ref. claim that such a negative branch in the heat capacity is a direct evidence of a first-order phase transition. Here we have a practical example that also a system with a second-order phase transition can show a negative specific heat and a backbending in the caloric curve, due to the non perfect equilibration. This effect increases with the size of the system and becomes stable in the limit $`N\mathrm{}`$ . This dynamical effect that simulates a first-order phase transition can be somehow explained by the presence of superdiffusion and Lévy walks in the out-of-equilibrium regime , that implies the coexistence of a liquid (clustered particles) and a gas (free particles) phase .

It is interesting to investigate the relaxation of the specific heat to the equilibrium value. In fig.4 we report $`C_V`$, calculated by means of eq. 4, versus time for the case U=0.6. $`C_V`$ is negative in the out-of-equilibrium regime, in correspondence of QSS, and converges very slowly to the equilibrium value. It is therefore necessary to wait for very long times, longer than for the temperature thermalization, to have also the thermalization of specific heat to positive values. A more complete study of the specific heat in HMF is in preparation. We notice in conclusions that the slow relaxation around the critical point persist notwithstanding the strong chaoticity found for the microscopic dynamics and reported in our previous work . The reason for that is not completely understood, though it has been found that the Lyapunov exponent is proportional to the kinetic fluctuations and some heuristic conjectures have been proposed .

## 4 CONCLUSIONS

We have presented equilibrium and out-of equilibrium microcanonical simulations of HMF, a simple system of interacting spins. Though the model is not explicitly built to describe nuclear system, we think it can give interesting information for the multifragmentation phase transition. In fact, HMF can be solved in the canonical ensemble and has a second-order phase transition in the thermodynamic limit. Moreover, the advantage with respect to Ising, lattice gas models and percolation, is that HMF allows a dynamical microcanonical approach to the study of phase transitions in a finite system. HMF has a very rich dynamical behavior in a transient out-of-equilibrium regime whose timescale depends on the energy and on the size of the system. When the microcanonical simulations are started in an out-of-equilibrium initial state, for example in a “water bag”, we find the appearance of quasi-stationary-states. In correspondence we have a caloric curve with a well defined backbending and negative specific heat. These indications of a first-order phase transition, at variance with the second-order phase transition predicted in the thermodynamic limit, are effects of the non-equilibrium. Certainly the non-equilibrium features of the model are strictly linked to the long-range nature of the interaction, however the general validity of these results is not completely clear and further work in this direction is needed.

We thank X. Campi, M. Pettini and C. Tsallis for useful discussions. Part of this work has been done in collaboration with S. Ruffo.