text:
# Generic Model of Morphological Changes in Growing Colonies of Fungi

## I Introduction

Since the pioneering work of D’Arcy Thompson  at the beginning of this century understanding the origin of complex structures involving living organisms has been a challenge for physicists and biologists alike. In particular the patterns formed in growing colonies of microorganisms, like bacteria  and fungi  have attracted much interest because they exhibit properties akin to those in inanimate systems far from equilibrium . Early mathematical models of fungal growth , considered fungi as purely additive assemblages of discrete individual hyphal units that duplicate at regular intervals and homogeneously. These models predict the crossover from exponential growth to linear growth of the peripheral zone as the colony matures. However, they are unable to account for many of the observable macroscopic properties (as morphological changes for instance) exhibited by fungi.

Fungal growth experiments on agar plates have shown  that filamental growth leads to a fractal colony with patterns that have much in common with those observed in diffusion-limited aggregation (DLA) processes . Similar branched morphologies have also been observed in bacterial growth  where the relevance of DLA was nicely demonstrated . Sophisticated communication mechanisms, including chemotactic signals and feedback response, were later proposed  in order to understand the origin of branched growth in bacteria. A similar mathematical modelling might be adequate to explain the DLA-like patterns of hyphae. However, quantitative measurements of the fractal dimension of mycelia  have shown that the dimension varies in time and from one species to another. Also, the growth mechanism of DLA clusters seems rather different from the hyphal extension of fungi. Some reaction-diffusion models have also been studied , which seem to capture esential features of the filamental fungi.

In contrast, under appropriate environmental conditions, colonies of many fungal species are also able to display a compact phase  with uni-cellular non-cooperative growth and generate non-DLA patterns. In this case, transitions between different morphologies are also observed  when the relevant parameters in the growth conditions are varied. A mathematical model able to include all these fungal patterns on solid media has remained elusive up to now . Here we show that a simple stochastic model containing some generic biological mechanisms is able to reproduce a variety of patterns observed in the growth of fungal colonies on solid agar plates. The model takes into account the effect of the accumulation of waste products, which play an important role as inhibitors for cell division.

Matsuura and Miyazima  carried out experiments with the fungus Aspergillus oryzae on agar plates of different stiffness and with different concentrations of nutrients (glucose). They observed that the fungal colonies on solid media formed a variety of shapes (see Figure 1). In nutrient-rich conditions, the colonies showed a homogeneous and smooth growth front (Fig. 1a). In contrast, growth in glucose-poor conditions resulted in colonies with rougher surfaces (Fig. 1b). The physiological ability of the fungus to divide and grow is diminished at lower temperatures , when even morphological instabilities (Fig. 1c) or an incipient branched growth (Fig. 1d) can be observed depending on the nutrient concentration. Similar patterns have recently been obtained in experiments with colonies of the yeast Pichia membranaefaciens on solidified agarose . How can these morphological changes be qualitatively described and understood within one simple mathematical model?.

## II The Model

It is known that metabolic products generated by fungal cells during the degradation of nutrients diffuse through the medium and may act as impurities inhibiting further cell division . We have considered the effect of toxic metabolites on the morphology of the colony. We have studied a model in which the growth probability at every site is coupled with an inhibitory field generated by the fungus itself. The model is defined on a two dimensional lattice as follows. At time $`t`$ a lattice site can be either occupied by a fungal cell or vacant. Two quantities are assigned to every site $`\stackrel{}{x}`$ in the lattice: (1) the total age of the nearest neighbour cells $`A(\stackrel{}{x},t)`$ of that site and (2) the field $`c(\stackrel{}{x},t)`$ containing the concentration of waste products. The spread of inhibitors from a single cell situated at $`\stackrel{}{x}_0`$ is modelled by means of a diffusion equation with a source:

$$\frac{c(\stackrel{}{x},t)}{t}=D^2c(\stackrel{}{x},t)+s\delta (\stackrel{}{x}\stackrel{}{x}_0)\mathrm{\Theta }(t\tau _0),$$
(1)
where $`\mathrm{\Theta }`$ is the step function ($`\mathrm{\Theta }(u)=1`$ for $`u0`$ and $`\mathrm{\Theta }(u)=0`$ for $`u<0`$). $`D`$ is the diffusion constant of the metabolites, $`s`$ is the inhibitor chemicals’ production rate, and $`\tau _0`$ is the time at which the polluting cell at $`\stackrel{}{x}_0`$ was born. The total density of waste products at any vacant site at time $`t`$ is calculated by adding the contributions coming from every occupied site in colony situated within a distance $`d`$.

Colony growth occurs because of the division of individual cells, thus only nearest neighbours of occupied sites have a chance of becoming occupied. We assume that the probability for a vacant site $`\stackrel{}{x}`$ of being occupied in the next time step increases with the total age $`A(\stackrel{}{x},t)`$ of the occupied nearest neighbours of that vacant site $`\stackrel{}{x}`$. Therefore, at every time step, we assign a growth probability $`P(\stackrel{}{x},t)`$ to every vacant site $`\stackrel{}{x}`$ given by

$$P(\stackrel{}{x},t)=F[\theta A(\stackrel{}{x},t)\varphi (c)],$$
(2)
which depends on the two fields $`A(\stackrel{}{x},t)`$ (total age of neighbouring occupied sites) and $`c(\stackrel{}{x},t)`$ (waste concentration) introduced above. In Eq.(2), $`F`$ is an arbitrary monotonous increasing function satisfying $`F(0)=0`$ and $`F(\mathrm{})=1`$ (properties required for $`P`$ to be a probability). The results do not depend on the detailed form of $`F`$ and in the following we use $`F(x)=\mathrm{tanh}(x)`$ for simplicity. The suppression of cell division due to the accumulation of metabolites is represented by $`\varphi (c)`$ that is a decreasing function of the concentration of waste chemicals. In order to see a significant effect of the inhibitors, the decay of $`\varphi `$ with $`c`$ must be rapid enough. In the simulations we are presenting in this paper $`\varphi (c)=\mathrm{exp}(c/c_0)`$ was assumed, $`c_0`$ being a concentration threshold above which the effect of inhibitors becomes important. Other simple functional forms may be used and similar patterns are obtained provided that $`\varphi (c)`$ decays rapidly enough. Finally, the external parameter $`\theta `$ controls directly the growth rate and is associated with the food supply (media richer in nutrients correspond to larger values of $`\theta `$).

## III Results

We carried out extensive simulations of the model. The simulations were performed on a two dimensional triangular lattice with open boundary conditions, simultaneously updating Eqs.(1) and (2) in every time step. The position of the front is determined and the height $`h(x,t)`$ at every column $`x`$ is used to study the front dynamics. In order to mimic the experiments, the growth is initiated from a straight line of randomly occupied sites at $`h(x,0)=0`$. The dynamics of the growing front is characterized by the height fluctuations, which are measured by means of the front width over the total system of size $`L`$, $`W(L,t)=\sqrt{(hh)^2}`$. The use of other quantities to measure the width is helpful when instabilities or overhangs appear, like for instance the difference $`Max[h]Min[h]`$ between the maximum and minimum heights, and gave similar results.

The model has two external parameters, namely: the growth rate $`\theta `$, and the concentration threshold $`c_0`$ for inhibitory effects. These are to be thought of as effective parameters and their relation to experimental variables such as temperature, concentration of nutrients, etc. may not be a simple one. It is to be noted that, if experimental conditions change, for instance in such a way that the absorption of inhibitors by the fungal cells is enhanced, the effect may be well represented by either a lower threshold $`c_0`$ or a wider spread distance $`d`$ (or even a combination of both).

We performed computer simulations of the model for different values of $`\theta `$ and $`c_0`$. In the following we discuss the different front morphologies that are obtained depending on the values of these two parameters, and how they are related to experimental observations. In Fig. 2 we present a summary of the resulting numerical fungal patterns.

### A High threshold for inhibitory effects

The first morphological transition occurs in the limit $`c_0\mathrm{}`$, in which the inhibitors become irrelevant ($`\varphi (c)1`$). We found that the surface of the colony is flat (Fig. 2a) when inhibitors have little effect (for $`c_01`$) and the medium is rich in nutrients. Conversely, the growth front becomes rougher and Eden-like (Fig. 2b) for a poorer medium. In fact, a quantitative measure of the roughness (see below) indicates that the fractal dimension of the front is $`1.5`$  as in the Eden model . These morphological phases are separated by a continuous phase transition at a critical value of $`\theta =\theta _c`$. For large growth rates (rich medium) $`\theta >\theta _c`$ the front is flat, $`W\mathrm{const}`$. In contrast, if the growth rate is reduced below $`\theta _c`$ (poor medium) the front becomes rougher and exhibits the dynamical scaling behaviour typical of a scale-free roughening process . As can be seen in Fig. 3, in the rough phase the width grows like a power law $`Wt^\beta `$ until a stationary regime is reached and $`WL^\alpha `$, where $`L`$ is the system size. Our determination of the time exponent $`\beta =0.24\pm 0.02`$ and the roughness exponent $`\alpha =0.46\pm 0.05`$ indicates that the model belongs to the Edwards-Wilkinson universality class  in this phase. The computation of the roughness exponent $`\alpha `$ allows us to determine the fractal dimension $`D_f`$ of the front through the relation $`D_f=2\alpha 1.54`$, which is consistent with Eden-like growth. This morphological transition from rough to flat growth has recently been studied in detail by us  for a much simpler version of this model, which corresponds to the limit $`c_0\mathrm{}`$. We found that $`\theta _c=0.183\pm 0.003`$ in the limit case of $`c_0=\mathrm{}`$. When approaching the roughening transition from the rough phase ($`\theta <\theta _c`$) to the flat phase ($`\theta >\theta _c`$), a diverging correlation lenght $`\xi |\theta \theta _c|^\nu `$ appears. The critical exponent was found to be $`\nu =1.10`$. At the phase transition the scaling exponents can be calculated by mapping the problem onto directed percolation (see Ref.  for further details). Experimental realizations of transitions from a flat to a rough front can been observed in the Matsuura and Miyazima experiments shown in Figure 1a-1b. As our model predicts the richer the medium the flater the front.

### B Low threshold for inhibitory effects

If the effect of toxic inhibitors on the fungus is enhanced, further interfacial instabilities in the front appear. Lower temperatures are less favorable for the physiological activity of the fungus  and its capacity to divide. One can model this effect by a smaller threshold $`c_0`$ for the inhibition of cell division. One can observe (Fig. 2c) that, as $`c_0`$ is decreased, deep grooves and overhangs develop in the colony and some regions with high concentration of inhibitors are eventually surrounded by the colony front. The presence of a unstable growth dynamics can be realized by the appearance of an extremely rapid (exponential) growth of the interface width $`W(L,t)`$ in this regime, as shown in Fig. 4. This regime is to be compared with the patterns observed in the expriments (Fig. 1c).

Finally, in the worst of the situations, the growth becomes very localized and only the dominant filaments branch out into the empty medium (Fig. 2d). This is due to the building up of extremely high concentrations of waste products in some regions that the fungus cannot invade any longer. There is qualitative agreement with the experiment shown in Fig. 1d.

The four major morphological phases are summarized in the qualitative phase diagram in Fig. 5. The solid line that separates the rough from the flat phase in Fig. 5 has recently been shown to be a nonequilibrium continuous phase transition and can be mapped into directed percolation . However, we were unable to find a sharp transition between the remaining phases.

## IV Discussion

Our results indicate that the morphologies generated during the growth of fungal colonies on solid agar can be modeled by taking into account (1) the building up of a concentration of waste inhibitors in the medium and (2) the provision of nutrients. We have studied a simple model which allows us to understand the interplay between these two mechanisms. The morphological tranformations that we have found can be compared with recent experiments on pattern formation in fungal growth.

In our model, a constant growth rate $`\theta `$ has been assumed. The effect of the exhaustion and diffusion of nutrients can be modeled by replacing $`\theta `$ by a diffusive field, which is depleted at a certain rate at sites occupied by fungal cells. This would have a major impact on the front dynamics as the front will eventually stop growing because of agar exhaustion.

Our model approach to pattern formation in fungal colony growth is stochastic and includes in a natural way the random nature of cell division. This together with the accumulation of inhibiting chemicals produces some of the different morphological phases observed in experiments. This is in contrast to other existing models based upon non-linear diffusion, in which evolution is deterministic and patterns arise from nonlinear instabilities. They represent a different level of description of the same phenomena. Both type of models coexist at different levels of description.

Finally, while our model can account for the morphological changes observed in some fungi growing in agar culture, we should mention that this model may not be appropiate for describing the morphological patterns observed in other fungal species. Our results compare well with the fungus Aspergillus oryzae which, like many members of this genus (and also the genus Penicillum), form dense colonies that expand only very slowly and seem to be self-inhibitory.

The authors are in debt to S. Matsuura for kindly providing us photographs of their experiments. We thank T. Sams for useful discussions on his work and M.J. Carlile, S. Matsuura, and E. Calle for helpful comments and correspondence. This work has been supported by a grant from the European Commission (Contract No. HPMF-CT-99-00133)).


format: text

text:
# Electroweak Knot

## Abstract

We demonstrate the existence of stable knot solitons in the standard electroweak theory whose topological quantum number $`\pi _3(S^2)`$ is fixed by the Chern-Simon index of the $`Z`$ boson. The electroweak knots are made of the helical magnetic flux tube of $`Z`$ boson which has a non-trivial dressing of the Higgs field, which could also be viewed as two quantized flux rings linked together whose linking number becomes the knot quantum number. We estimate the mass of the lightest knot to be around $`21TeV`$.

electroweak knot, topological knot

Ever since Dirac proposed his theory of monopoles the topological objects in physics have been the subject of intensive studies dirac ; skyr . In particular the finite energy topological solitons have been widely studied in theoretical physics skyr ; thooft . A remarkable type of solitons is the knots which have recently appeared almost everywhere, in nuclear physics in Skyrme theory fadd1 ; cho01 ; cho1 , in plasma physics in coronal loops fadd2 , in condensed matter physics in one-gap as well as multi-gap superconductors cho2 ; cho3 ; baba , and in atomic physics in two-component Bose-Einstein condensates cho4 ; ruo .

Of these, the Faddeev-Niemi knot in Skyrme theory and the superconducting knot in ordinary superconductor are particularly important. The Faddeev-Niemi knot is the prototype of the all the knots, which comes from the Skyrme-Faddeev Lagrangian

$`_{SF}={\displaystyle \frac{\mu ^2}{2}}(_\mu \widehat{n})^2{\displaystyle \frac{1}{4}}(_\mu \widehat{n}\times _\nu \widehat{n})^2.`$ (1)
The Lagrangian has the following equation cho01

$`\widehat{n}\times ^2\widehat{n}{\displaystyle \frac{g}{\mu ^2}}(_\mu H_{\mu \nu })_\nu \widehat{n}=0,`$
$`H_{\mu \nu }={\displaystyle \frac{1}{g}}\widehat{n}(_\mu \widehat{n}\times _\nu \widehat{n}),`$ (2)
which admits not only the knot but also the helical baby skyrmion, a twisted baby skyrmion which is periodic in $`z`$-coordinate. The importance of the helical baby skyrmion is that the Faddeev-Niemi knot is nothing but a vortex ring made of the helical baby skyrmion with two periodic ends smoothly connected together. This tells that the Faddeev-Niemi knot originates from the helical baby skyrmion cho1 . The knot is non-Abelian, because the knot topology $`\pi _3(S^2)`$ comes from the $`SU(2)`$ symmetry of the theory.

On the other hand, the superconducting knot in ordinary superconductor is an Abelian knot whose knot topology is fixed by the Chern-Simon index of the electromagnetic potential cho2 . Nevertheless this knot is closely related to the other knots in a fundamental way. It can also be viewed as a vortex ring made of an helical magnetic vortex, the twisted Abrikosov vortex cho2 . This tells that the existence of a helical vortex is an essential condition for a knot. It guarantees the existence of the knot.

The purpose of this Letter is to demonstrate the existence of an electroweak knot in Weinberg-Salam theory. We show that the electroweak knot is made of two quantized neutral magnetic flux rings linked together, the first one winding the second $`m`$ times and the second one winding the first $`n`$ times, whose linking number $`mn`$ becomes the knot quantum number. Furthermore we predict that the lightest knot has mass around $`21TeV`$ and size of $`3.5\times 10^{18}m`$. We also show that the knot has both topological and dynamical stability, which comes from the twisted topology of the neutral magnetic field. If confirmed by experiments, the knot could constitute the first topological particle in high energy physics.

To establish the existence of the electroweak knot we first need to understand the deep connection which exists between the Skyrme theory and the non-Abelian gauge theory. Let ($`\widehat{n}_1,\widehat{n}_2,\widehat{n}`$) be a right-handed orthonormal basis in $`SU(2)`$ space, and consider the following decomposition of the potential $`\stackrel{}{A}_\mu `$ into the restricted potential $`\widehat{A}_\mu `$ and the valence potential $`\stackrel{}{X}_\mu `$ cho80 ; cho81 ,

$`\stackrel{}{A}_\mu =A_\mu \widehat{n}{\displaystyle \frac{1}{g}}\widehat{n}\times _\mu \widehat{n}+\stackrel{}{X}_\mu =\widehat{A}_\mu +\stackrel{}{X}_\mu ,`$
$`A_\mu =\widehat{n}\stackrel{}{A}_\mu ,\stackrel{}{X}_\mu =X_\mu ^1\widehat{n}_1+X_\mu ^2\widehat{n}_2.`$ (3)
Notice that $`\widehat{A}_\mu `$ is precisely the connection which leaves $`\widehat{n}`$ invariant under parallel transport,

$`\widehat{D}_\mu \widehat{n}=_\mu \widehat{n}+g\widehat{A}_\mu \times \widehat{n}=0.`$ (4)
Under the infinitesimal gauge transformation

$`\delta \widehat{n}=\stackrel{}{\alpha }\times \widehat{n},\delta \stackrel{}{A}_\mu ={\displaystyle \frac{1}{g}}D_\mu \stackrel{}{\alpha },`$ (5)
one has

$`\delta A_\mu ={\displaystyle \frac{1}{g}}\widehat{n}_\mu \stackrel{}{\alpha },\delta \widehat{A}_\mu ={\displaystyle \frac{1}{g}}\widehat{D}_\mu \stackrel{}{\alpha },`$
$`\delta \stackrel{}{X}_\mu =\stackrel{}{\alpha }\times \stackrel{}{X}_\mu .`$ (6)
This tells that $`\widehat{A}_\mu `$ by itself describes an $`SU(2)`$ connection which enjoys the full gauge degrees of freedom. Furthermore $`\stackrel{}{X}_\mu `$ forms a gauge covariant vector field under the gauge transformation. But what is really remarkable is that the decomposition is gauge-independent. Once $`\widehat{n}`$ is given, the decomposition follows automatically independent of the choice of a gauge cho80 ; cho81 .

Notice that $`\widehat{A}_\mu `$ retains the full topological characteristics of the original non-Abelian potential. Clearly, $`\widehat{n}`$ defines $`\pi _2(S^2)`$ which describes the non-Abelian monopole cho80m . Besides, with the $`S^3`$ compactification of $`R^3`$, $`\widehat{n}`$ describes the Hopf invariant $`\pi _3(S^2)\pi _3(S^3)`$ which characterizes both topologically distinct vacua and instanton number cho79 ; cho02 . Furthermore $`\widehat{A}_\mu `$ has a dual structure,

$`\widehat{F}_{\mu \nu }=(F_{\mu \nu }+H_{\mu \nu })\widehat{n}\text{,}`$
$`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu \text{,}`$
$`H_{\mu \nu }={\displaystyle \frac{1}{g}}\widehat{n}(_\mu \widehat{n}\times _\nu \widehat{n})=_\mu C_\nu _\nu C_\mu ,`$ (7)
where $`A_\mu `$ and $`C_\mu `$ are the (chromo)electric and (chromo)magnetic potential cho80 ; cho81 . Notice that $`H_{\mu \nu }`$ here is exactly the same two-form appeared in (2), which admits the potential $`C_\mu `$ because it is closed.

The decomposition (3) reveals the deep connection between the Yang-Mills theory and the Skyrme theory. To see this notice that with

$`\stackrel{}{H}_{\mu \nu }=_\mu \stackrel{}{C}_\nu _\nu \stackrel{}{C}_\mu +g\stackrel{}{C}_\mu \times \stackrel{}{C}_\nu =H_{\mu \nu }\widehat{n},`$
$`\stackrel{}{C}_\mu ={\displaystyle \frac{1}{g}}\widehat{n}\times _\mu \widehat{n},`$ (8)
we have cho02

$`_{SF}={\displaystyle \frac{1}{4}}\stackrel{}{H}_{\mu \nu }^2{\displaystyle \frac{\mu ^2}{2}}\stackrel{}{C}_\mu ^2.`$ (9)
This tells that the Skyrme-Faddeev theory can be interpreted as a massive Yang-Mills theory where the gauge potential has the special form (8). This is a first indication that the Weinberg-Salam theory could admit a knot similar to Faddeev-Niemi knot. Our decomposition (3), which has recently been referred to as the “Cho decomposition” fadd3 ; gies ; zucc , plays a crucial role in QCD, in particular in the calculation of the effective action of QCD cho02c .

With these preliminaries we now demonstrate the existence of an electroweak knot. Consider the Weinberg-Salam Lagrangian

$`={\displaystyle \frac{1}{4}}\stackrel{}{F}_{\mu \nu }^2{\displaystyle \frac{1}{4}}G_{\mu \nu }^2|\stackrel{~}{D}_\mu \varphi |^2+m^2\varphi ^{}\varphi {\displaystyle \frac{\lambda }{2}}(\varphi ^{}\varphi )^2,`$
$`\stackrel{~}{D}_\mu \varphi =(_\mu +{\displaystyle \frac{g}{2i}}\stackrel{}{\sigma }\stackrel{}{A}_\mu +{\displaystyle \frac{g^{}}{2i}}B_\mu )\varphi .`$ (10)
With the decomposition (3) we can identify the Higgs field $`\rho `$ and $`W`$ boson $`W_\mu `$ by,

$`\varphi ={\displaystyle \frac{\rho }{\sqrt{2}}}\xi (\xi ^{}\xi =1),W_\mu ={\displaystyle \frac{1}{\sqrt{2}}}(X_\mu ^1+iX_\mu ^2),`$ (11)
and express (10) in terms of the physical fields alone

$`={\displaystyle \frac{1}{2}}(_\mu \rho )^2{\displaystyle \frac{g^2}{4}}\rho ^2W_\mu ^{}W_\mu {\displaystyle \frac{g^2+g^2}{8}}\rho ^2Z_\mu ^2`$
$`{\displaystyle \frac{\lambda }{8}}\left(\rho ^2\rho _0^2\right)^2{\displaystyle \frac{1}{4}}F_{\mu \nu }^{(\mathrm{em})}{}_{}{}^{2}{\displaystyle \frac{1}{4}}Z_{\mu \nu }^2`$
$`{\displaystyle \frac{1}{2}}|(D_\mu ^{(\mathrm{em})}W_\nu D_\nu ^{(\mathrm{em})}W_\mu )+ie{\displaystyle \frac{g}{g^{}}}(Z_\mu W_\nu Z_\nu W_\mu )|^2`$
$`+ieF_{\mu \nu }^{(\mathrm{em})}W_\mu ^{}W_\nu +ie{\displaystyle \frac{g}{g^{}}}Z_{\mu \nu }W_\mu ^{}W_\nu `$
$`+{\displaystyle \frac{g^2}{4}}(W_\mu ^{}W_\nu W_\nu ^{}W_\mu )^2,`$ (12)
where $`\rho _0^2=2m^2/\lambda `$, $`D_\mu ^{(\mathrm{em})}=_\mu +ieA_\mu ^{(\mathrm{em})}`$, $`A_\mu ^{(\mathrm{em})}`$ and $`Z_\mu `$ are the electromagnetic potential and the $`Z`$ boson. We emphasize that this expression is obtained without any gauge fixing, which is made possible with the gauge independent decomposition (3).

Now, with the ansatz

$`A_\mu ^{(\mathrm{em})}=0,W_\mu =0,`$ (13)
we have the following equation

$`^2\rho g_0^2Z_\mu ^2\rho ={\displaystyle \frac{\lambda }{2}}(\rho ^2\rho _0^2)\rho ,`$
$`_\mu Z_{\mu \nu }=j_\nu =g_0^2\rho ^2Z_\nu ,`$ (14)
where $`g_0=\sqrt{g^2+g^2}/2`$ and $`Z_{\mu \nu }=_\mu Z_\nu _\nu Z_\mu `$. This admits not only a helical vortex but also a knot made of a twisted magnetic field of the $`Z`$ boson.

We construct the helical vortex first. Choose the cylindrical coodinates $`(\varrho ,\phi ,z)`$ and the following ansatz

$`\rho =\rho (\varrho ),`$
$`Z_\mu ={\displaystyle \frac{1}{g_0}}\left(nZ_1(\varrho )_\mu \phi +mkZ_2(\varrho )_\mu z\right).`$ (15)
With this we have

$`Z_{\varrho \phi }={\displaystyle \frac{n}{g_0}}\dot{Z}_1,Z_{\varrho z}={\displaystyle \frac{mk}{g_0}}\dot{Z}_2,Z_{\phi z}=0,`$
$`j_\mu =g_0\rho ^2\left(nZ_1_\mu \phi +mkZ_2_\mu z\right),`$ (16)
which clearly shows that both the magnetic flux and the supercurrent have a helical structure.

With the ansatz (14) is reduced to

$`\ddot{\rho }+{\displaystyle \frac{1}{\varrho }}\dot{\rho }\left({\displaystyle \frac{n^2}{\varrho ^2}}Z_1^2+m^2k^2Z_2^2\right)\rho ={\displaystyle \frac{\lambda }{2}}(\rho ^2\rho _0^2)\rho ,`$
$`\ddot{Z}_1{\displaystyle \frac{1}{\varrho }}\dot{Z}_1g_0^2\rho ^2Z_1=0,`$
$`\ddot{Z}_2+{\displaystyle \frac{1}{\varrho }}\dot{Z}_2g_0^2\rho ^2Z_2=0.`$ (17)
Now, we impose the following boundary condition

$`\rho (0)=0,\rho (\mathrm{})=\rho _0,Z_1(0)=1,Z_1(\mathrm{})=0.`$ (18)
As for $`Z_2`$ we require that the vortex carries a non-vanishing supercurrent. This uniquely fixes the boundary condition for $`Z_2`$, which dictates that $`Z_2(\mathrm{})=0`$ but $`Z_2(0)`$ must have a logarithmic divergence.

With this we obtain the helical vortex shown in Fig.1. Notice that, when $`m=0`$, the solution (with $`Z_2=0`$) describes the well-known $`Z`$ boson vortex vach . But when $`m`$ is not zero, it describes a helical vortex which has a non-vanishing supercurrent (not only around the $`z`$-axis but also) along the $`z`$-axis,

$`i_{\widehat{z}}=mkg_0{\displaystyle \frac{}{}}{\displaystyle \rho ^2Z_2\varrho 𝑑\varrho 𝑑\phi }={\displaystyle \frac{2\pi mk}{g_0}}(\varrho \dot{Z}_2)|_{\varrho =0}.`$ (19)
This confirms that the logarithmic divergence of $`Z_2`$ at the origin is what we need to make the vortex a superconducting string. But here, of course, the supercurrent is neutral, not electromagnetic.

Clearly the vortex has the magnetic field $`H_{\widehat{z}}`$ and the quantized magnetic flux $`\mathrm{\Phi }_{\widehat{z}}`$ along the $`z`$-axis

$`H_{\widehat{z}}={\displaystyle \frac{n}{g_0}}{\displaystyle \frac{\dot{Z}_1}{\varrho }},\mathrm{\Phi }_{\widehat{z}}={\displaystyle \frac{}{}}{\displaystyle H_{\widehat{z}}\varrho 𝑑\varrho 𝑑\phi }={\displaystyle \frac{2\pi n}{g_0}}.`$ (20)
But it also has the magnetic field $`H_{\widehat{\phi }}`$ around the vortex

$`H_{\widehat{\phi }}={\displaystyle \frac{mk}{g_0}}\dot{Z}_2.`$ (21)
Unfortunately this produces an infinite the magnetic flux $`\mathrm{\Phi }_{\widehat{\phi }}`$ around the vortex, because of the singularity at the origin. So one might conclude that the helical vortex is unphysical which does not exist, because one need an infinite energy to create it.

The importance of the helical vortex, however, is not in that it is physical but in that it ensures the existence of a knot. To see this notice that with the vortex we can make a vortex ring by smoothly bending and connecting two periodic ends. In this vortex ring the infinite magnetic flux $`Z_2`$ becomes finite because the finite supercurrent (19), in the vortex ring, should produce a finite magnetic flux passing through the disk of the ring. This tells that the infinite magnetic flux $`\mathrm{\Phi }_{\widehat{\phi }}`$ of the vortex is an artifact of a straight vortex, which disappears in the vortex ring. Furthermore, we can certainly make the finite magnetic flux passing through the knot disk to have the quantized value $`2\pi m/g_0`$ by adjusting the supercurrent of the ring with $`k`$. Remarkably, this vortex ring now becomes a topologically stable knot. To see this notice that the magnetic flux of the knot is made of two flux rings, a $`2\pi m/g_0`$ flux around the knot tube and a $`2\pi n/g_0`$ flux along the knot. Moreover the two flux rings are linked together, whose linking number becomes $`mn`$. This is precisely the mathematical description of a knot, two rings linked together. This assures that the vortex ring does indeed become a topological knot. The knot quantum number is described by the Chern-Simon index of $`Z_\mu `$ field,

$`Q={\displaystyle \frac{g_0^2}{32\pi ^2}}{\displaystyle ϵ_{ijk}Z_iZ_{jk}d^3x}=mn,`$ (22)
which describes the non-trivial topology $`\pi _3(S^2)`$ of the magnetic field $`Z_{\mu \nu }`$. Notice that this is formally identical to the quantum number of the Faddeev-Niemi knot cho01 ; cho1 . The only difference is that here the (chromo)magnetic potential $`C_\mu `$ is replaced by $`Z_\mu `$. Furthermore, just as in the Faddeev-Niemi case, the Chern-Simon index is given by the linking number of two magnetic fluxes.

Obviously the two flux rings linked together can not be separated with a continuous deformation of the field configuration. This provides the topological stability of the knot.

Furthermore, this topological stability of the knot is backed up by a dynamical stability. This is because the supercurrent along the knot now generates a net angular momentum around the knot which naturally provides the centrifugal repulsive force to prevent the collapse of the knot. This tells that the knot is dynamically stable. Another way to understand the dynamical stability is to notice that, when the knot tries to shrink, the energy density of magnetic field trapped in the knot disk increases inevitably. This generates a repulsive force against the collapse, which makes the knot stable.

We emphasize hat the stability of the knot crucially depends on the helical structure of the supercurrent and magnetic field. Without this there is neither the topological stability nor the dymamical stability.

We can estimate the mass of the electroweak knot. To do this notice that in the absence of the Higgs field the knot becomes almost identical to the Faddeev-Niemi knot in Skyrme theory. From this observation we may estimate the energy density (per length) and the energy, thus the radius of the knot as cho4 ; ussr

$`=\pi n\rho _0^2,E16\pi ^2\times 3^{3/8}(mn)^{3/4}g_0\rho _0,`$
$`R={\displaystyle \frac{E}{2\pi }}{\displaystyle \frac{8g_0\times 3^{3/8}}{\rho _0}}{\displaystyle \frac{m^{3/4}}{n^{1/4}}}.`$ (23)
So, with the experimental values $`g_00.36`$ and $`\rho _0246GeV`$ we expect the lightest electroweak knot to have mass around $`21TeV`$ and radius of about $`3.5\times 10^{18}m`$, with tube size $`0.7\times 10^{18}m`$.

Clearly the electroweak knot is closely related to the Faddeev-Niemi knot. But one can not overemphasize the striking similarity between the electroweak knot and superconducting knot in ordinary superconductor cho2 . Mathematically they are identical. The only difference is that the superconducting knot exists at the atomic scale and carries the real electric supercurrent, whereas the knot here exists at the electroweak scale and carries the neutral current. It is really remarkable that mathematically identical knots can exist in totally different physical enviornments, in $`eV`$ scale and in $`TeV`$ scale.

We believe that our analysis has established the existence of a electroweak knot beyond reasonable doubt. Of course, one might like to see an analytic solution of the knot. Unfortunately this is impossible. At present even the simplest knot does not allow an analytic solution fadd1 . But one can construct an approximate solution analytically which has all the charactristic features of the electroweak knot, using the toroidal coordinates and adopting an educated ansatz cho2 ; cho5 . This provides another evidence which endorses the fact that the electroweak knot is for real. A challenging task now would be to confirm the existence of the electroweak knot by high energy experiments. We hope that the LHC at CERN could confirm the existence of the electroweak knot.

The physical implications of the electroweak knot and the details of our analysis will be published separately cho5 .

ACKNOWLEDGEMENT

We thank Professor C. N. Yang for the illuminating discussions. The work is supported by the Basic Research Program of Korea Science and Engineering Foundation (Grant R02-2003-000-10043-0) and by the BK21 project of the Ministry of Education.


text:
# Spectrum Results with Kogut-Susskind Quarks

## Acknowledgements

I am grateful to the organizers of Lattice-01 for the opportunity to present this talk. I thank Claude Bernard, Joachim Hein, Takashi Kaneko, Francesco Knechtli, Kostas Orginos, and Howard Trottier for providing results used in this talk. This work was supported by the U.S. DOE.


text:
# Brightest Cluster Galaxy Formation in Low-Redshift Galaxy Clusters

## 1. Introduction

Galaxy clusters provide a laboratory in which one can study the effects of galaxy formation and evolution over time. A detailed understanding of low-redshift clusters would furnish us with a benchmark in which to compare observations of higher redshift galaxy clusters. Cluster properties, such as luminosity functions, have the potential to yield valuable clues concerning the dynamical evolution of cluster galaxies, including the formation of Brightest Cluster Galaxies (BCGs).

Recent studies have concentrated on the merging process as an efficient way of creating BCGs (Merritt 1984; Schombert 1988). The formation of cD-like galaxies (we define cD galaxies as those with an extended halo (e.g., Schombert 1988) as well as those galaxies which are halo-dominated; Brown 1997) has been problematic since the merging process cannot account for the presence of an extended halo (Merritt 1984; Tremaine 1990).

From a survey of 45 low-redshift ($`0.04z0.18`$) galaxy clusters, López-Cruz et al. 1997 proposed that the envelope of cD galaxies was created through the process of tidal disruption of dwarf galaxies (e.g., dwarf spheroidals). Evidence for this hypothesis was obtained by examining cluster luminosity functions (LFs). López-Cruz et al. (1997) found that clusters which contain a cD galaxy (Bautz Morgan types I, I-II; Rood-Sastry type “cd”) had flatter faint-end slopes (as parameterized by the Schechter function; Schechter 1976) than those which did not. Further support for this model was found by examining the dwarf-to-giant galaxy ratio (D/G), which was found to increase with increasing Bautz Morgan type (López-Cruz 1997).

## 2. Cluster Survey

In order to test the dwarf galaxy disruption model, we have conducted an imaging survey of low-redshift ($`0.02z0.04`$) Abell clusters using the KPNO 0.9m telescope + 8k mosaic camera. We have obtain R band images for 27 clusters, of which a sub-sample of 11 were also observed in the B band. The 8k mosaic camera has a field-of-view of approximately one square degree and allows us to cover a region 2 to 4 Mpc in length, centred on each cluster.

The purpose of this survey is to test the dwarf galaxy disruption hypothesis by; 1) sampling cluster luminosity functions to a depth adequate to measure the faint-end slope, 2) measuring cluster LFs over a wide area in order to sample the LF as a function of cluster-centric radius, 3) determining the colour gradient of the BCG envelope and compare it to the radial colour gradient of the dwarf galaxy population, and 4) measuring the properties (LFs, spatial distributions, etc.) of morphologically selected galaxy sub-populations.

In this paper we concentrate on results based on the analysis of cluster luminosity functions (we use $`H_o=50kms^1Mpc^1`$ and $`q_o=0`$ throughout).

### 2.1. Data Reductions and Photometry

Object detection, classification, and photometry was performed using PPP (Yee 1991). All bright, early-type galaxies were modeled and subtracted from their parent image using the same procedure as in Brown (1997). This allowed us to accurately measure the magnitudes of overlapping fainter galaxies which becomes a serious problem in the cluster core. Complete photometric catalogs of each cluster contain approximately 25,000 objects.

Since we do not have redshifts for every cluster galaxy, luminosity functions are created from the statistical subtraction of background galaxies. The background correction was determined by observing six controls fields, covering a total of six square degrees in area.

## 3. Luminosity Functions

The direct comparison of LFs for various surveys has been problematic. Many studies have relied on photographic plates and have compared cluster LFs measured using a different radial cutoff from the BCG (e.g., Lugger 1986). Since we expect that environment may have an influence on the overall galaxy population mixture (morphology-density relation; Dressler 1980), it is prudent to compare cluster LFs using the same dynamical radius. For this study, we have chosen to measure our LFs using the $`r_{200}`$ radius (the radius within which the average density is 200 times the critical density). Since we lack radial velocity measurements, we use the correlation between $`r_{200}`$ and Bgc (a richness indicator; Yee & López-Cruz 1999) as measured for the CNOC 1 cluster survey (Carlberg, Yee, & Ellingson 1997). We find that as cluster richness increases, the value of $`r_{200}`$ also increases. Using this relationship, we can estimate a value of $`r_{200}`$ for a given value of Bgc (which we can calculate from our photometry). In this way, we are able to make a more direct and robust comparison of cluster LFs.

In an effort to improve the signal-to-noise of cluster galaxies over fore/background galaxies, we impose a colour criteria from which we select cluster members. Cluster colour-magnitude diagrams clearly show the well known colour-magnitude relation (CMR; Baum 1959) for early-type cluster galaxies. The colour-magnitude diagram can reveal the presence of background contaminating clusters (López-Cruz 1997). We have used the colour-magnitude diagram to establish an upper envelope in which to reject redder background galaxies.

### 3.1. Composite Luminosity Function

We computed the composite LF for two radial bins centred on the BCG. The inner area was chosen to have a radius of $`r_{200}`$, and the outer bin had an inner radius of $`r_{200}`$ and an outer radius $`2r_{200}`$. Clusters whose photometry was 100% complete to $`M_R=16.0`$, and which had complete imaging out to a radius of $`2r_{200}`$, were included in the final sample. This sample consisted of 17 clusters, 7 of which have colour information which allowed us to impose a colour selection criteria based on the cluster’s CMR. All magnitudes were K corrected based on Coleman, Wu, & Weedman (1980) and corrected for galactic extinction (Burstein & Heilis 1982).

Individual cluster LFs were constructed through the statistical subtraction of background galaxy counts as measured from six control fields. The rejection criteria used for each cluster was also applied to the background counts (including normalization with respect to cluster area). Each cluster LF was then summed to produce a composite LF for each radial bin.

Figure 1 shows the composite LF computed for the inner and outer radial regions. It is clear from this figure that a single Schechter function is inadequate to describe the composite LFs. We have therefore used the sum of two Schechter functions to fit the LFs in both cases. For the inner radial bin ($`rr_{200}`$), we find a faint-end slope of $`\alpha _{inner}=1.81\pm 0.02`$. The faint-end slope for the outer radial bin ($`r_{200}r2r_{200}`$) was measured to be $`\alpha _{outer}=2.07\pm 0.02`$. It is clear from this result that the inner region has a significantly flatter faint-end slope than the outer region.

## 4. Discussion and Conclusions

Schechter (1976) described an empirical function which was shown to provide an adequate fit to a sample of 14 cluster LFs. Many studies have investigated the possibility that the LF is “universal” (e.g., Lugger 1986). The results of several recent studies (e.g., López-Cruz et al. 1997; Driver, Couch, & Phillipps 1998), clearly demonstrate that cluster LFs are not universal but are instead a function of environment (this is really just a re-statement of the density-morphology relation).

The comparison of the composite LFs in this study support the conclusion that cluster LFs are not universal and are in fact dependent upon cluster-centric radius (ie., environment). This result is also consistent with the dwarf galaxy disruption model since one expects that tidal forces acting on dwarf galaxies would be greater near the central regions of clusters. A reduction in the fraction of dwarf galaxies, relative to the remaining galaxy cluster population, would be expected for the inner cluster regions.

## References

Baum, W.A. 1959, PASP, 71, 106

Brown, J.P. 1997, Ph.D. thesis, University of Toronto

Burstein, D., & Heiles, C. 1982, AJ, 87, 1165

Carlberg, R.G., Yee, H.K.C., & Ellingson, E. 1997, ApJ, 478, 462

Coleman, G.D., Wu, Chi-Chao, & Weedman, D.W. 1980, ApJS, 43, 393

Dressler, A. 1980, ApJ, 236, 351

Driver, S.P., Couch, W.J., & Phillipps, S. 1998, MNRAS, 301, 369

López-Cruz, O. 1997, Ph.D. thesis, University of Toronto

López-Cruz, O., Yee, H.K.C., Brown, J.P., Jones, C., & Forman, W. 1997, ApJ, 476, L97

Lugger, P.M. 1986, ApJ, 303, 535

Merritt, D. 1984, ApJ, 276, 26

Schechter, P.L. 1976, ApJ, 203, 297

Schombert, J.M. 1988, ApJ, 328, 475

Tremaine, S.D. 1990, in Dynamics and Interactions of Galaxies, Wielen, R., Ed. (Heidelberg FRG: Springer-Verlag), 394

Yee, H.K.C. 1991, PASP, 103, 396

Yee, H.K.C., & López-Cruz, O. 1999, AJ, 117, 1985


text:
# The overlap operator as a continued fractionBased on a talk by U. Wenger

## 1 Introduction

Let us start with the explicit form of the massive overlap Dirac operator with mass $`\mu `$,

$$D(\mu )=\frac{1}{2}(1+\mu )+\frac{1}{2}(1\mu )\gamma _5\text{sgn}(H(m)),$$
(1)
where $`H(m)`$ denotes some hermitian Dirac operator with mass $`m<0`$ of the order of the inverse lattice spacing. One method to compute the action of the sign-function on a vector is by approximating the function by a rational polynomial ,

$$\text{sgn}(x)\text{sgn}_{\{n,m\}}(x)\frac{P_n(x)}{Q_m(x)},$$
(2)
where $`P_n`$ and $`Q_m`$ are irreducible polynomials of degree $`n`$ and $`m`$, respectively. The reason for using rational approximations is that for a given approximation error in the infinity-norm the required degree for rational approximations is usually much smaller than the one for polynomial ones, although both converge exponentially in their degree.

There are various methods to perform such an approximation. The methods differ in how they interpret the notion of the goodness of the approximation. One can use any definition of a norm characterising the deviation of two functions from each other.

The infinity-norm, for example, is defined as the maximal error over a given interval $`[a,b]`$. There exists an algorithm (Remez algorithm ) to find the best (rational) approximation $`r(x)`$ to a function $`f(x)`$ with respect to the infinity-norm (see also the work of Zolotarjov ). That is, the approximation $`r(x)`$ to a given function $`f(x)`$ is the rational function of a given degree that minimises the maximum value of the quantity $`|f(x)r(x)|`$ over the interval $`[a,b]`$.

In order to approximate the sign-function by a rational polynomial we choose to find an approximation for a constant odd function over the interval $`[\epsilon ,1]`$. To do so we approximate $`\frac{1}{\sqrt{x}}`$ on the interval $`[\epsilon ^2,1]`$ to obtain the rational approximation $`r(x)`$ and consider the function $`xr(x^2)`$.

The resulting approximation can then be written as

$$\text{sgn}_{\{2n+1,2m\}}(x)=xc_0\frac{_{k=1}^n(x^2+p_k)}{_{k=1}^m(x^2+q_k)}$$
(3)
with real and positive $`p_k`$ and $`q_k`$. This factorised form is most stable for numerical evaluation.

One advantage of using rational polynomials is that one can use the partial fraction decomposition, e.g. for $`n=m`$,

$$\text{sgn}_{\{2n+1,2n\}}(x)=x\left(c_0+\underset{k=1}{\overset{n}{}}\frac{c_k}{x^2+q_k}\right)$$
(4)
which allows us to compute the action of the sign-function at roughly the cost of a single conjugate gradient (CG) inversion of $`H^2`$ .

## 2 Continued fractions

It is well known that rational functions are closely related to finite (truncated) continued fractions which in turn play an important role in the analytic theory of rational functions . Indeed, we can write any rational approximation to the sign-function as a continued fraction

$$\text{sgn}_{\{2n+1,2n\}}(x)=\beta _0x+\frac{\alpha _1}{\beta _1x+{\displaystyle \frac{\mathrm{}}{\mathrm{}+{\displaystyle \frac{\alpha _{2n}}{\beta _{2n}x}}}}}.$$
(5)
Then we observe that the inverse of eq. (5) is just the $`(1,1)`$-component of the Schur complement of a larger system

$$\left(\begin{array}{ccccc}\beta _0x& \sqrt{\alpha _1}& & & \\ \sqrt{\alpha _1}& \beta _1x& & & \\ & & \mathrm{}& & \\ & & & \beta _{2n1}x& \sqrt{\alpha _{2n}}\\ & & & \sqrt{\alpha _{2n}}& \beta _{2n}x\end{array}\right).$$
(6)
Instead of inverting (4) using a two-level nested CG(-like) procedure we can now find the inverse of (6), i.e. its (1,1)-component, in one go, using just one CG(-like) search in a single higher dimensional Krylov space .

The operation naturally associated with rational functions is the multiplication and inversion which, of course, is most trivial. For the overlap operator, eq. (1), matters are not so simple since the additional term proportional to $`\gamma _5`$ does not commute with $`\text{sgn}(H)`$ and its approximations. We can, however, still make use of some properties of (finite) continued fractions. The most important one is contained in the following theorem on equivalence transformations: Two continued fractions $`\{\alpha _0;\{\alpha _1,\beta _1\},\mathrm{},\{\alpha _n,\beta _n\}\}`$ and $`\{\alpha _0^{};\{\alpha _1^{},\beta _1^{}\},\mathrm{},\{\alpha _n^{},\beta _n^{}\}\}`$ are equivalent if and only if there exists a sequence of non-zero constants $`c_n`$ with $`c_0=1`$ such that

$`\alpha _n^{}=`$ $`c_nc_{n1}\alpha _n,`$ $`n=1,2,3,\mathrm{},`$
$`\beta _n^{}=`$ $`c_n\beta _n,`$ $`n=0,1,2,\mathrm{}.`$
This can most easily be seen by writing down the continued fraction after a transformation,

$$\beta _0+\frac{c_1\alpha _1}{c_1\beta _1+{\displaystyle \frac{c_1c_2\alpha _2}{c_2\beta _2+{\displaystyle \frac{\mathrm{}}{\mathrm{}+{\displaystyle \frac{c_{n1}c_n\alpha _n}{c_n\beta _n}}}}}}}.$$
(7)
In fact the coefficients $`c_i`$ parametrise the equivalence class of a continued fraction corresponding to a given rational polynomial. It is important to note that although the equivalence transformation leaves the value of the continued fraction invariant it affects the spectrum and correspondingly the condition number of the larger system (6). One might therefore hope to find within an equivalence class of continued fractions the one with the smallest condition number.

## 3 Application to the overlap operator

It is easy to apply the equivalence transformation to an overlap operator in which the sign-function is approximated by a rational function and expanded as a continued fraction. In complete analogy to the previous section we express the linear system of the hermitian overlap Dirac operator $`\gamma _5D(\mu )\psi =b`$ as a simple five dimensional system $`H^{(5)}(\mu )\mathrm{\Psi }=\chi `$ where $`\mathrm{\Psi }=(\psi ,\varphi _1,\mathrm{},\varphi _{2n+1})^T`$, $`\chi =(b,0,\mathrm{},0)^T`$ are five dimensional fermion fields and $`H^{(5)}(\mu )`$ is a block tridiagonal matrix similar to (6) with $`A_+\gamma _5,c_1k_1H,+c_2k_2H,\mathrm{},c_{2n+1}k_{2n+1}H`$ on the diagonal, $`\sqrt{c_1A_{}},\sqrt{c_1c_2},\mathrm{},\sqrt{c_{2n}c_{2n+1}}`$ with trivial Dirac and colour structure on the off diagonals and $`A_\pm =1/2(1\pm \mu )`$. The $`k_i`$’s are uniquely determined by the rational approximation of the sign-function while the $`c_i`$’s are free parameters describing the equivalence class of continued fractions associated with the given rational approximation.

We observe that an equivalence transformation on the continued fraction acts in a similar way to block Jacobi preconditioning. Indeed we can write $`H^{(5)}(\mu )=CH^{(5)}(\mu ;c_i=1)C`$ where $`C`$ is a block diagonal matrix having the four dimensional unit matrix $`\mathrm{𝟏}`$ in the first block and $`\sqrt{c_i}\mathrm{𝟏},i=1,\mathrm{},2n+1`$ in the remaining ones. For the linear system under consideration this preconditioning comes for free since the first component of $`\mathrm{\Psi }`$ and $`\chi `$ are not affected at all.

It turns out that within an equivalence class the condition number of $`H^{(5)}`$, $`\kappa (H^{(5)})`$, can vary by orders of magnitude. An obvious choice for the equivalence transformation parameters is $`c_i=1/k_i`$, and indeed, with this choice the matrix is already well-conditioned. In order to find the optimal values for the $`c_i`$ in the sense of giving the lowest $`\kappa (H^{(5)})`$ one can study the matrix in the free case where the spectrum can be calculated analytically. A non-linear minimisation of $`\kappa (H^{(5)})`$ then yields a set of optimal values for the coefficients $`c_i`$.

How much is the effective gain in a practical application? In order to address this question we generated an ensemble of pure gauge field configurations with the Wilson gauge action at $`\beta =5.9`$ on $`6^4`$ lattices. We used the standard Wilson Dirac operator $`H_W(m)`$ with $`m=1.5`$ as the input operator for the overlap and a $`n=8`$ rational approximation for the inverse square root on the interval $`[\epsilon ^2,1]=[0.00724,1]`$. This generates a five dimensional operator with an extent of 18 in the fifth dimension. Then we calculated the inverse of the hermitian overlap operator $`\gamma _5D(\mu )`$ using the five and four dimensional formulation with different quark masses $`\mu `$ on a random vector.

Figure 1 shows the average number of $`H_W`$-matrix multiplications for the inversion in three cases: the pole approximation, eq. (4), ($`4d`$ overlap), the five dimensional operators $`H^{(5)}`$ with $`c_i=1/k_i`$ ($`5d`$ matrix) and $`H^{(5)}`$ with the optimal $`c_i`$’s as obtained from the free field approximation ($`5d`$ matrix optimised).

The results indicate that it is indeed possible to generate five dimensional systems which are well conditioned and that the obstacle of having large condition numbers for the five dimensional system can be avoided. For small quark masses the inversion seems to behave better than the one of the corresponding four dimensional system. This is of particular importance in the context of dynamical overlap fermion simulations where the inversion of the fermion matrix plays a crucial role . Another advantage of the five dimensional formulation using continued fractions is its obvious suitability for parallel computations.


text:
# Lower limits on the maximum orbital frequency around rotating strange stars

## 1 Introduction

The recently discovered kHz quasi periodic oscillations may be used as a probe of the inner regions of accretion disks in compact stars and hence also of some of the properties of the central object. For example, these oscillations have been used to derive estimates of the mass of the central source in some sources (e.g., Kaaret et al., 1997; Zhang et al., 1998; Kluźniak 1998). All these estimates assume that the maximum observed QPO frequency is the orbital frequency at the innermost stable circular orbit (ISCO), as suggested by Kluźniak et al. (1990).

In the first study of the possible presence of strange stars in those LMXBs which exhibit kHz QPOs, Bulik et al. (1999a,b) showed that slowly rotating strange stars, described by the simple MIT bag model with massless and non-interacting quarks, have orbital frequencies at the marginally stable orbit higher than the maximum frequency of 1.07 kHz in 4U 1820-30, reported by Zhang et al. (1998). However, the ISCO frequencies can be as low as 1 kHz when more sophisticated models of quark matter (MIT bag with massive strange quarks and lowest order QCD interactions) and/or rapid stellar rotation are taken into account (Stergioulas et al., 1999; Zdunik et al., 2000a,b). The main result of these considerations was that the lowest orbital frequency at the ISCO is attained either for non-rotating massive configurations close to their maximum mass limit, or for configurations at the Keplerian mass-shedding limit. In order to obtain an orbital frequency in the ISCO as low as 1.07 kHz for slowly rotating strange stars one has to consider a specific set of (rather extreme) MIT bag model parameters (Zdunik et al. 2000a, b) allowing gravitational masses higher than 2.2 $`M_{}`$. In contrast, for stars rotating close to the Keplerian limit, a low orbital frequency at the ISCO can be obtained for a broad range of stellar masses (Stergioulas et al., 1999).

A new description of quark matter (Dey et al., 1998) has been used for constructing numerical models of static strange stars and strange stars rotating with two different fixed frequencies 200 Hz and 580 Hz, and for calculating frequencies of marginally stable orbits around these stars (Datta et al., 2000). The authors conclude that very high QPO frequencies in the range of $`1.9`$ to $`3.1`$ kHz would imply the existence of a non-magnetized strange star in an X-ray binary rather than a neutron star. Gondek-Rosińska et al. (2001) calculated the maximum orbital frequencies for the Dey et al. quark stars rotating at various rotation rates and showed that the rotational and maximum orbital frequencies for these stars are indeed much higher than those for neutron star models and for strange stars in the MIT bag model. The maximum orbital frequencies for the compact objects modeled with the Dey et al. (1998) equation of state are $`1.5`$ kHz, and always higher than the kHz QPO frequencies observed to date.

Recently, new complications in the problem were pointed out. The innermost stable circular orbit exists also in Newtonian gravity (Kluźniak, Bulik, Gondek-Rosińska 2000; Amsterdamski et al., 2000; Zdunik & Gourgoulhon 2001). The presence of a crust has been investigated for normal evolutionary sequences by Zdunik, Haensel & Gourgoulhon (2001). Typically the crust increases the maximum orbital frequency at the Keplerian limit.

In this paper we consider constant baryon mass sequences of strange stars described by the simplest MIT bag model (with massless and non-interacting particles) for various values of the baryon mass and calculate the frequency in the co-rotating ISCO (if it exists) for each sequence. We find and discuss lower limits to the maximum orbital frequency (taking into account various secular and dynamical stability limits, as well as physical constrains relevant for accretion in LMXBs). We then extend our results to all MIT equations of state. Our calculations are described in section 2 and the results are presented in section 3. Finally, section 4 contains a summary of our results.

## 2 The model star calculations

Typically, strange stars are modeled (Alcock et al., 1986; Haensel et al., 1986) with an equation of state based on the MIT-bag model of quark matter, in which quark confinement is described by an energy term proportional to the volume (Farhi and Jaffe 1984). The equation of state is given by the simple formula

$$P=a(\rho \rho _0)c^2,$$
$`(1)`$
where $`P`$ is the pressure, $`\rho `$ the mass-energy density, and $`c`$ the speed of light. It was shown by Zdunik (2000) that all MIT bag models can be approximated very well by equation (1), where $`a`$ is a constant of model-dependent value in the range 0.289 (for a strange quark mass of $`\mathrm{m}_\mathrm{s}`$= 250 MeV) to 1/3 (for massless quarks). In the numerical calculations reported in the present paper we describe strange quark matter using the simplest MIT bag model (with massless and non-interacting quarks), accordingly we take $`a=1/3`$ and $`\rho _0=4.2785\times 10^{14}\mathrm{g}\mathrm{cm}^3\rho _{\mathrm{fid}}`$ (corresponding to a bag constant of $`B=60\mathrm{MeV}\mathrm{fm}^3`$). In general, eq. (1) corresponds to self-bound matter with density $`\rho _0`$ at zero pressure and with a fixed sound velocity ($`\sqrt{a}c`$) at all pressures. For a fixed value of $`a`$, all stellar parameters are subject to scaling relations with appropriate powers of $`\rho _0`$ (see, e.g., Witten 1984; Zdunik 2000). Approximate scalings with $`a`$ were considered by Stergioulas et al. (1999) for the limiting case of rotation at Keplerian frequencies. Although a precise scaling-law with $`a`$ for other cases has not been derived, we will show that the lower limits on the maximum ISCO orbital frequency obtained for the simplest MIT bag model can easily be used to find lower limits on the maximum orbital frequency for any other value of $`a`$ and $`\rho _0`$.

We have calculated the maximum orbital frequency around uniformly rotating strange stars using two different highly accurate relativistic codes. The first code is based on the multi-domain spectral methods developed by Bonazzola et al. (1998). This method has been used previously for calculating properties of rapidly rotating strange stars described by the MIT bag model (Gourgoulhon et al., 1999) and by the Dey model (Gondek-Rosińska et al., 2000a). This code was also used for calculating the maximum orbital frequency around normal sequences of MIT strange stars (Zdunik et al., 2000b, Zdunik & Gourgoulhon 2001) and of the Dey model strange stars (Gondek-Rosińska et al., 2001). The second code (Stergioulas & Friedman 1995, see Stergioulas 1998 for a description) has been used to compute accurate constant baryon mass sequences and models at the Keplerian mass-shedding limit, where the first code has difficulties converging to a solution. In this code the equilibrium models are obtained following the KEH method (Komatsu et al., 1989), in which the field equations are converted to integral equations using appropriate Green’s functions. Models computed by both codes agree very well in all computed properties.

We calculate the maximum frequency of a co-rotating test particle in a stable circular orbit in the equatorial plane about a strange star, $`f_{\mathrm{orb},\mathrm{max}}`$. By testing stability of the orbital motion, we determine whether or not stable orbits extend to the surface of the star (at radius $`R_{}`$). If they do, the maximum orbital frequency is the keplerian frequency at the surface, $`f_{\mathrm{orb},\mathrm{max}}=f_{\mathrm{orb}}(R_{})`$. If they do not, we find the radius of the innermost marginally stable orbit, $`R_{\mathrm{ms}}`$, and its frequency $`f_{\mathrm{ms}}`$ (see, e.g. Friedman, Ipser & Parker 1986). In this latter case, the orbital frequency in the ISCO is the maximum circular orbital frequency, $`f_{\mathrm{orb},\mathrm{max}}=f_{\mathrm{ms}}`$ . This frequency can then be compared with the observed kHz QPOs, having in mind some model of QPOs relating their frequency to the orbital frequency and the location of the ISCO (e.g., Kluźniak et al., 1990, Abramowicz and Kluźniak 2001, Wagoner et al., 2001). In so doing, we neglect the effects of fluid pressure, magnetic fields, radiation drag, etc., on the location of the innermost circular orbit in the accretion disk.

We construct equilibrium sequences of rotating compact strange stars with constant baryon mass. We compute both normal and supramassive sequences. A sequence is called normal if it terminates at the zero angular momentum limit with a static, spherically symmetric solution, and it is called supramassive if it terminates at the axisymmetric stability limit, instead. The boundary between these two sequences is the sequence with the maximum baryon mass of a static configuration. The angular momentum (and the central density) of a star changes monotonically along each sequence.

For a supramassive constant baryon mass sequence a model is stable to axisymmetric perturbations, if $`\left({\displaystyle \frac{J}{\rho _\mathrm{c}}}\right)_{M_\mathrm{b}}<0`$ (or $`\left({\displaystyle \frac{M}{\rho _\mathrm{c}}}\right)_{M_\mathrm{b}}<0`$) where $`J`$ is the angular momentum, $`M`$ is the gravitational mass, $`M_\mathrm{b}`$ is the baryon mass of the star and $`\rho _\mathrm{c}`$ is the central density of the star (Friedman, Ipser & Sorkin 1988; Cook, Shapiro & Teukolsky 1994).

The star in a sequence reaches the Keplerian limit (the mass-shedding limit) if the velocity at the equator of a rotating star is equal to the velocity of an orbiting particle at the surface of the star.

## 3 Results

The properties of the ISCO around strange stars differ from those around neutron stars (Miller, Lamb & Cook 1998)—for strange (i.e., quark) stars of any mass the ISCO always exists at sufficiently rapid stellar rotation rate (Stergioulas et al., 1999; Zdunik et al., 2000b), and for strange stars of moderate and high baryon masses the ISCO always exists for any rotation rate (Zdunik et al., 2000b; Gondek-Rosińska et al., 2001). For lower mass stars (e.g., the $`1.6M_{}`$ baryon mass sequence in Fig. 1) stable orbits extend down to the stellar surface for moderate rotation rates, but the gap is present for both slowly and rapidly rotating stars. For lower masses yet, the gap is present only at high rotation rates and the ISCO is present even if relativistic effects are negligible (as in substellar-mass configurations, e.g., for $`M0.01M_{}`$). In the Newtonian limit, the gap between the marginally stable orbit and the stellar surface is produced by the oblateness of the rapidly rotating low-mass quark star (Kluźniak, Bulik & Gondek-Rosińska, 2000; Amsterdamski et al., 2000; Zdunik & Gourgoulhon 2001). For Keplerian configurations the ISCO is always well above the stellar surface for all possible masses.

We present the dependence of the maximum orbital frequency on the frequency of rotation $`f_{\mathrm{rot}}=\mathrm{\Omega }/2\pi `$ for a set of constant baryon mass sequences of strange stars in Figure 1. In Figure 2 we show how the gravitational mass depends on rotational frequency for the considered sequences. The normal sequences start at $`f_{\mathrm{rot}}=0`$, and continue with increasing angular momentum up to the Keplerian configuration. As it turns out, the very low mass limit, represented in Figure 1 by the case of $`M_\mathrm{b}=0.01M_{}`$, is described well by a uniform-density fluid, i.e., the Maclaurin spheroid (Chandrasekhar 1969). For this sequence, an ISCO exists only for rotation rates larger than $`970\mathrm{Hz}`$, and we indicate the onset of the (Newtonian) dynamical instability to non-axisymmetric perturbations by an asterisk.

The intermediate-mass case is shown by the constant baryon mass sequence with the baryon mass $`M_\mathrm{b}=1.6M_{}`$, for which the ISCO does not exist between the rotation rates with frequency of $`f_{\mathrm{rot}}300`$ Hz and $`f_{\mathrm{rot}}800`$ Hz (see also Zdunik et al., 2000b).

More massive stars in normal sequences are shown for the cases of $`M_\mathrm{b}=2.0M_{}`$ and $`2.3M_{}`$. In these cases, the marginally stable orbit exists along the entire constant baryon mass sequence. The ISCO frequency initially increases with rotational frequency, as calculated in the slow-rotation approximation by Kluźniak & Wagoner 1985. When the small-rotation-rate assumption is no longer satisfied, the ISCO frequency reaches a maximum and begins to decrease mainly due to the rotational deformation of the star (Sibgatulin & Sunyaev 1998; Shibata & Sasaki 1999). The sequence terminates when the star has reached the Keplerian configuration (mass-shedding limit). An extreme case of normal sequences is the sequence starting with the maximum mass static configuration of $`M_\mathrm{b}=2.624M_{}`$ (indicated as 2.6 $`M_{}`$ in Figures 1 and 2).

Supramassive sequences start at the onset of the instability to axisymmetric perturbations (thick dashed line in Figure 1 and in subsequent figures.) Sequences with $`M_\mathrm{b}=2.7M_{}`$ and $`2.8M_{}`$ are examples of low mass supramassive sequences. The behavior of the maximum orbital frequency vs. frequency of rotation is basically similar to the normal sequences described above.

Finally, we present the sequences of models with $`M_\mathrm{b}=3.2,\mathrm{\hspace{0.17em}3.4},`$ and $`3.6M_{}`$. For these models of extremely high mass, the rotational frequency decreases with increasing angular momentum because the stars are initially highly deformed. The ISCO frequency only decreases along these sequences until the stars reach the mass shedding limit. The configuration with the highest rotation frequency is not Keplerian. For a high-mass supramassive strange star sequence (e.g., $`M_\mathrm{b}=3.2`$), the Keplerian configuration is the one with the lowest rotational frequency in the sequence (see Fig. 2).

Note that the rotational frequency does not necessarily change monotonically along a constant baryon mass sequence, since a star changes its shape and, consequently, its moment of inertia with increasing angular momentum (see also Cook, Shapiro & Teukolsky 1992; Zdunik et al., 2000b; Gondek-Rosińska et al., 2000a). For a lower-mass supramassive sequence and a normal sequence the rotational frequency decreases with increasing angular momentum just below the mass-shedding limit (characteristic turning back in Figure 1 and 2 for each sequence). A similar behavior is observed in the case of models with small mass, e.g. for the constant baryon mass sequence with $`M_\mathrm{b}=0.01M_{}`$. Thus, the effect of spin-down with increasing angular momentum is of Newtonian origin, and is related to the change of shape of the star.

In Figure 1 one can see that the value of the ISCO frequency alone is not a sufficient indicator for the mass of the compact object. For a rapidly rotating strange star the same value of orbital frequency at the innermost stable circular orbit, e.g. 1.2 kHz, can correspond to a range of stars differing in mass by a factor of a hundred, or even more, as stressed already by Kluźniak, Bulik & Gondek-Rosińska 2000. However, it is worth noting that for rotational frequencies exceeding about 1.3 kHz, i.e., for rapidly rotating supramassive stars, the rotational period alone constrains the stellar mass to a narrow range.

### 3.1 Lower limits on the ISCO frequency

Figure 1 shows that at any given stellar rotational frequency the maximum orbital frequency $`f_{\mathrm{orb},\mathrm{max}}(f)`$ is limited from below by the least frequency among the following three “extreme” sequences of stars: a) the marginally stable configuration to axisymmetric perturbations (thick dashed line); b) configurations rotating at Keplerian frequencies (thick solid line) and c) the Newtonian limit of low-mass stars, tending to Maclaurin spheroids when $`M0`$, $`J0`$, $`\rho \rho _0`$ (solid line ending with the asterisk). For all these limiting cases, an ISCO exists. The maximum orbital frequencies are also bounded from above by a fourth limiting line (not shown), which will be discussed elsewhere: the envelope of highest values of $`f_{\mathrm{orb},\mathrm{max}}`$.

a) We find that at a given low or moderate value of stellar rotational frequency the lowest ISCO frequency is obtained for the configurations at the axisymmetric instability limit. As was shown by Gourgoulhon et al. (1999) and Gondek-Rosińska et al. (2000a) a supramassive strange star, just as a supramassive neutron star (Cook et al., 1992), prior to its collapse to a black hole increases its spin as it loses angular momentum. In such cases, two models with the same baryon mass but with different central density and angular momentum can have the same rotational frequency (see inset in Figure 1). Along a sequence of fixed rotational frequency, the model with the lowest ISCO frequency is not the model with the highest gravitational or highest baryon mass. The dependence of the ISCO frequency as function of the stellar central density for sequences of stars with a fixed rotational frequency is shown in Figure 3. As the central density increases, the ISCO frequency decreases and the sequence passes through the configuration with maximum gravitational mass (denoted with short dashed line in Figure 3 and Figure 2), and then the maximum baryon mass configuration (denoted with the dash-dotted line). Before reaching the minimal value of $`f_{\mathrm{orb},\mathrm{max}}`$, the sequence reaches the limit of stability to axisymmetric perturbations (thick dashed line). Thus, for a fixed rotational frequency, the lowest possible maximum orbital frequency is attained for the model which is marginally stable to gravitational collapse.

b) At high rotation rates the lowest ISCO frequency is attained for the Keplerian (mass-shedding) configurations (thick solid line in Figure 1). The stellar configurations at this limit are highly deformed. We calculated Keplerian configurations beginning with the fastest rotating stable configurations (with $`M_\mathrm{b}=3.76M_{}`$), ending with a strange star with $`M_\mathrm{b}=0.4M_{}`$. Note that for rapidly rotating (close to the Keplerian limit) strange stars, the ISCO frequency is smaller than the rotational frequency.

c) In the small-mass limit the ISCO only exists at high stellar rotation rates. This was discovered numerically, and the Newtonian origin of the ISCO in this case was first understood by comparison with toy models of disk-like distributions of matter (Kluźniak, W., Bulik, T.,& Gondek-Rosińska 2000; Zdunik & Gourgoulhon 2001). In the Newtonian limit ($`M0`$, $`J/M^2\mathrm{}`$) strange stars tend to a nearly uniform-density configuration ($`\rho \rho _0`$), and these can be described well by the Maclaurin spheroids, for which analytic formulae for orbital frequencies have been derived and found to agree with the numerical models (Amsterdamski et al., 2000). The functional dependence of the maximum orbital frequency on the stellar rotation rate is nearly universal for all strange stars with masses lower than 0.1 $`M_{}`$, in the model considered here.

### 3.2 The lower limits on ISCO for other MIT bag models

Figure 1 and the discussion above assumed the simplified MIT bag model with massless and non-interacting quarks. The results can easily be scaled for any value of the density at zero pressure, with the use of the scaling $`f\rho _0^{1/2}`$, where $`f`$ denotes either of the frequencies (rotational or orbital), and $`M\rho _0^{1/2}`$. We can easily find lower limits on ISCO for any other MIT bag model since Zdunik (2000) showed that all MIT bag model equations of state can be well approximated by eq. (1).

In order to find the dependence of the lower limits on $`a`$ we display in Figure 4 the limits calculated for $`\mathrm{m}_\mathrm{s}=0`$ and for $`\mathrm{m}_\mathrm{s}=250`$MeV. We have chosen these values of the strange quark mass because it corresponds to the highest value $`a=0.333`$, and the lowest value of $`a=0.289`$ for the MIT bag model equation of state. The dependence on the density of quark matter at zero pressure has been scaled out by multiplying the frequencies with the square root of $`\rho _{\mathrm{fid}}/\rho _0`$, where our fiducial value of $`\rho _0`$ was used: $`\rho _{\mathrm{fid}}=4.2785\times 10^{14}\mathrm{g}\mathrm{cm}^3`$. The numerical values of the frequencies in Fig. 4 correspond to the value of the bag constant $`B=60\mathrm{MeV}\mathrm{fm}^3`$ used throughout this paper.

The limit on the maximum orbital frequency, set by the low mass stars does not depend on the value of $`a`$. The limit at rapid rotation rates depends only slightly on the value of $`a`$. The main difference appears in the slow and moderate rotation regions, where the lower limit is set by the sequence of marginally stable stars. However, here one can use the fact that the limit is only shifted along the vertical axis. We find that the formula $`f_{\mathrm{orb},\mathrm{max}}(a,f_{\mathrm{rot}})=0.3f_{\mathrm{rot}}+f_{\mathrm{orb},\mathrm{max}}(a,0)`$ is a good approximation for describing the lower limit on the ISCO frequency for slowly and moderately rotating strange stars described by any MIT equation of state. Here, the maximum orbital frequency of the nonrotating model is given by its Schwarzschild-metric value

$$f_{\mathrm{orb},\mathrm{max}}(a,0)=2.2\mathrm{kHz}/(M(a)/\mathrm{M}_{})$$

(e.g., Kluźniak et al., 1990), and $`M(a)`$ is the maximum mass in general relativity of a static configuration of the (quark) fluid described by eq. (1)—in short, the maximum mass of a non-rotating strange star.

### 3.3 Factors that could alter the lower limits

There are several factors which could affect the lower limits considered in the previous section. We discuss here: i) different kinds of instabilities that could operate in rapidly rotating strange stars, ii) the existence of a solid crust on the accreting strange stars, iii) details of the accretion-induced spin-up of strange stars in LMXBs.

i) Rapidly rotating compact stars may be subject to non-axisymmetric rotational instabilities. A widely used indicator for the onset of instability is the ratio of the kinetic energy to the absolute value of the gravitational potential energy $`T/|W|`$. We plot the value of $`T/|W|`$ (defined as in Friedman et al., 1986) for several constant baryon mass sequences in Figure 5. The values of $`T/|W|`$ for strange stars rotating at the Keplerian velocity are very large and can be even higher than $`T/|W|=0.32`$.<sup>1</sup><sup>1</sup>1In Gourgoulhon et. al., 1999 (see their Fig. 7), $`T/|W|<0.24`$, because configurations denoted as Keplerian, were in fact configurations at which the numerical procedure stopped converging. The $`T/|W|`$ is very high for all linear self-bound equations of state (see Gondek-Rosińska et al., 2000b for different MIT bag models and Gondek-Rosińska et al., 2000a for the strange stars described by the Dey model). The values of $`T/|W|`$ for strange stars are much larger than for any neutron star models (for which the maximum $`T/|W|`$ ranges typically from 0.1 to 0.14, (see e.g. Nozawa et al., 1998). In contrast to neutron stars, the $`T/|W|`$ for strange stars increases with decreasing stellar mass for all models. Very large values of $`T/|W|`$ are characteristic not only for supramassive stars but also for rapidly rotating strange stars with small and moderate masses. The large value of $`T/|W|`$ results from a flat density profile combined with strong equatorial flattening of rapidly rotating strange stars. No stability analysis in general relativity has been performed yet for strange stars.

There exists two different classes and corresponding time-scales for non-axisymmetric instabilities: a dynamical instability, growing on the hydrodynamical time-scale and secular instabilities (due to dissipative mechanisms) which grow on much longer time-scales. In Newtonian theory, rapidly rotating Maclaurin spheroids are secularly unstable to a bar-mode formation, if $`T/|W|>(T/|W|)_{\mathrm{crit}}=0.1375`$ (Jacobi/Dedekind bifurcation point) and dynamically unstable if ($`T/|W|>0.2738`$). A secular instability can grow only in the presence of dissipative mechanisms, such as shear viscosity (Roberts & Stewartson 1963) or gravitational radiation (the Chandrasekhar-Friedman-Schutz instability (hereafter CFS instability; Chandrasekhar 1970; Friedman & Schutz 1978; Friedman 1978). If both viscosity and gravitational radiation are operative, they act against each other. Which of the two dissipative mechanisms dominates is highly dependent on the temperature of the star. By a suitable choice of the strength of viscosity relative to gravitational radiation, it is possible to stabilize the Maclaurin sequence up to the point of dynamical instability (Shapiro & Teukolsky 1983).

In general, the viscosity-driven instability could be operating in accreting (old and cold) compact objects, while the CFS instability could operate in both newly-born compact stars and in LMXBs. The most important pulsation modes for the CFS instability are the $`r`$\- and $`f`$-modes (see Stergioulas 1998, Andersson & Kokkotas 2001, Friedman and Lockitch 2001, for recent reviews of these instabilities).

The dynamical instability against bar-mode deformation of rapidly rotating polytropic stars was studied by Shibata, Baumgarte & Shapiro (2000) in full general relativity and by Saijo et al., 2001 in the first post-Newtonian approximation. The authors considered differentionally rotating proto-neutron stars. Their main result is that in general relativity the onset of dynamical instability occurs for a somewhat smaller critical value of $`(T/|W|)_{\mathrm{crit}}0.240.25`$ than the Newtonian value $`0.27`$ for incompressible Maclaurin spheroids.

The few numerical studies of secular instabilities in rapidly rotating stars, in general relativity show that nonaxisymmetric modes driven unstable by viscosity no longer coincide with those driven unstable by gravitational radiation. The $`(T/|W|)_{crit}`$ depends not only on the dissipative mechanism but also on the compaction parameter $`M/R`$ of a star. The viscosity-driven instability for neutron stars has been studied in general relativity by Bonazzola, Frieben & Gourgoulhon (1996, 1998), who find that it sets in for $`T/|W|`$ somewhat larger than the Newtonian value of $`(T/|W|)_{\mathrm{crit}}=0.1375`$. A stronger weakening of the viscosity-driven instability by general relativity is suggested by a post-Newtonian study of the viscosity-driven instability in incompressible stars by Shapiro & Zane (1998), who find that the $`(T/|W|)_{\mathrm{crit}}`$ could be as large as 0.25 for $`M/R0.2.`$ On the other hand, the gravitational-radiation driven instability is strengthened by relativity and the bar ($`l=2`$) $`f`$-mode can become unstable for typical neutron star equations of state (Stergioulas & Friedman 1998). For a wide range of realistic equations of state, Morsink, Stergioulas & Blattning (1999) find that the critical $`T/|W|`$ for the onset of the instability in the $`l=2`$ $`f`$-mode is only $`T/|W|0.08`$ for models with mass $`1.4M_{}`$ and $`T/|W|0.06`$ for maximum mass models.

The CFS-instability of $`r`$-modes in strange stars has been studied by Madsen (2000), in the Newtonian limit. According to this study, the instability can operate at a minimum rotational frequency of $`400`$Hz and at a temperature of $`T10^7`$K. Since strange stars can, at most, have only a thin solid crust, they cannot store heat in the interior during accretion. Thus, accreting strange stars are likely to have interior temperatures comparable to $`T10^7`$K. If the $`r`$-mode instability is not damped by other mechanisms then the above results seem to rule out the presence of low and moderate mass stars with a normal quark matter in LMXBs (if the observed 1.07 kHz QPO is the orbital frequency at the ISCO).

If, contrary to our current understanding, the $`r`$-mode instability is not effective in accreting strange stars, then the $`f`$-mode instability or the viscosity-driven instability could be more relevant. More detailed studies of these instabilities specifically for relativistic strange stars are needed. The dynamical instability, even if present at $`T/|W|0.24`$, only marginally affects the limits presented in Fig. 1.

ii) Another factor affecting the lower limits presented in Figure 1 would be the presence of a solid crust. This effect has been discussed by Zdunik, Haensel & Gourgoulhon (2001). The presence of the crust will not affect the lower limit in the small and moderate rotation case. Here the stability of the strange stars is basically determined by the stability of the strange core and the properties of the crust are negligible. However, for rapid rotation, when the lower limit on the maximum orbital frequency is determined by the Keplerian configurations, the presence of the crust will increase the lower limit. The crust is of low density and for rapidly rotating strange stars it becomes extended in the equatorial region. Thus the mass-shedding limit is reached at lower rotation rates than in the case of bare strange stars. Zdunik, Haensel & Gourgoulhon (2001) showed that the values of maximum orbital frequency for strange stars with maximum crust (up to the neutron drip point), rotating at the mass-shedding limit, are about $`20\%`$ larger compared to case of bare strange stars. Also, for strange stars with a crust, the ISCO exists only for models with large mass.

iii) When a strange star is spun-up by accretion of material that leaves the ISCO and plunges onto the surface of the star, the largest attainable rotation rate is not precisely the rotation rate at the mass-shedding limit. Material that leaves from the ISCO arrives at the surface of a strange star after following a path in which it conserves angular momentum. At the surface, the material has angular velocity larger than what it had at the ISCO, but still somewhat smaller than the angular velocity needed to spin-up the star to the angular velocity of the equilibrium configuration at the mass-shedding limit. The difference between the maximum rotation and mass-shedding configurations is, however, small.

## 4 Summary

Relativistic orbital frequencies are crucially related to the spin-up of neutron stars and black holes and to high energy accretion phenomena in low mass X-ray binaries, including quasi-periodic oscillations in the X-ray flux. We compute the maximum orbital frequency of stable circular motion around strange stars, and present it as a function of the stellar rotation rate for normal and supramassive constant baryon mass sequences of stars rotating at all possible rates. We find that in contrast to neutron stars, the mass of a rapidly rotating strange star cannot be even approximately inferred from the orbital frequencies alone. For strange stars described by the MIT bag model, the same frequency in the innermost stable orbit (e.g., 1.2 kHz) is obtained for strange stars with rotational periods ranging from infinity to about 0.6 ms, and the mass ranging from that of a planetoid to about three solar masses. One should note however that the rapidly rotating stars in this region are quite likely to be unstable.

We find that for slow to moderate rotation rates the lower limit on the maximum orbital frequency around strange stars is attained for configurations marginally stable to axisymmetric perturbations, while at rapid rotation rates it is reached for configurations at the mass-shedding limit. There are several factors which could affect the lower limit for rapidly rotating strange stars, such as different instabilities or the existence of a solid crust. The $`r`$-mode instability has the potential to rule out strange stars with small and moderate masses in LMXBs (assuming that the observed 1.07 kHz QPO in 4U 1820-30 is close to the Keplerian frequency at the innermost stable orbit of an accretion disk around a strange star). If the $`r`$-mode instability is not effective, then the $`f`$-mode instability and the viscosity-driven instability may become relevant. On the other hand, the dynamical instability only marginally affects our computed lower limits on the maximum orbital frequency. Our results apply to all MIT bag models.

###### Acknowledgements.

This work has been funded by the following grants: KBN grants 5P03D01721, 2P03D00418, 2P03D01816; the Greek-Polish Joint Research and Technology Programme EPAN-M.43/2013555 and the EU Programme ”Improving the Human Research Potential and the Socio-Economic Knowledge Base” (Research Training Network Contract HPRN-CT-2000-00137).


text:
# The BTeV Vertex Trigger System1footnote 11footnote 1Presented at the Snowmass Summer Study on the Future of Particle Physics, Snowmass, Colorado, June 30–July 21, 2001.

## I Introduction

BTeV is an approved Fermilab experiment that will use the Tevatron Collider to study heavy-quark physics. BTeV data-taking is planned to commence in Tevatron Run IIB (starting about 2006). To achieve the highest possible sensitivity to mixing, CP violation, and rare decays, a very ambitious trigger system is envisioned that will examine every beam crossing for evidence of secondary vertices.

The trigger is conceived in three levels. The Level 1 trigger (the lowest trigger level) will employ custom pattern-recognition hardware to perform fast trackfinding on every beam crossing, using information from the silicon-pixel vertex detector BTeV-proposal ; FPIX2 . Events for which the pixel detectors show evidence of secondary vertices will be passed on to Level 2, which will use information from additional tracking detectors to improve the vertex resolution and refine the Level 1 decision. Events passing Level 2 will undergo a more complete event reconstruction at Level 3, with information available from all detectors. (In addition we are building a high-$`p_t`$-muon trigger, which we do not discuss further here; similar hardware will be used to implement both triggers.)

## II Vertex trigger strategy

BTeV is designed for a 7.6-MHz beam crossing rate at an average of 2 interactions per crossing. In reconstructing tracks in real time from all crossings, the key problem is the combinatoric processing time that typically limits the speed of trackfinding algorithms. Since the initial conception of BTeV BTeV-EoI , our thinking has evolved on the best approach for meeting this challenge. A key part of the solution is to optimize the design of the vertex detector not only for good vertex resolution, but also for fast pattern recognition. For this reason we are designing thin silicon pixel detectors with fast, data-driven readout FPIX2  and locating them within the spectrometer magnet, so that momentum information is available already at the lowest trigger level (see Fig. 1). Fast readout with high hit-rate capability requires rectangular pixels (Fig. 2), to maintain excellent measurement resolution in at least one coordinate while providing sufficient room behind each pixel for the complex pixel signal-processing and readout circuitry.

Our initial conception of the trigger algorithm Kaplan-B99 , based on work by the Penn group Husby , required that each pixel station by itself determine a short track segment or “mini-vector.” The mini-vector should have sufficient resolution to select only a small search region in the next station, with mean occupancy well below 1 hit per crossing. To accomplish this, we 1) provided each station with three pixel planes and 2) divided each plane logically into 32 pie-shaped sectors (“$`\varphi `$-slices”). Monte Carlo simulation gave a mean hit occupancy of 0.23/interaction/$`\varphi `$-slice. By providing in each station two “precision-$`y`$” planes and one “precision-$`x`$” plane, we could determine, using the hits in a single station, a mini-vector accurate to $`\pm 1`$mr in $`y`$ and $`\pm 50`$mr in $`x`$ that would thus project into a very small fraction ($``$ $`2\times 10^5`$) of the area of the next station, within which the occupancy due to other tracks or random hits would be negligibly small.

This vertex-detector strategy led to a high-performance, cost-effective trigger implementation BTeV-pTDR ; Kaplan-B99 . However, it required an extra pixel plane per station compared to the minimum needed for off-line event reconstruction. In view of the total thickness of the triplet-scheme vertex detector (approaching 1 radiation length), as well as its cost, a 2-plane-per-station solution was clearly desirable, and in 1999 we embarked on a two-year exploration of alternative trigger approaches that would be compatible with 2 planes per station. Two-plane stations required giving up the mini-vector approach; the key trigger problem then becomes the trackfinding combinatorics.

We have studied two alternative approaches compatible with two-plane stations. At 2 interactions/crossing, a pixel-plane quadrant has 8 hits on average; the probability of $`>`$16 hits is $`<`$1%. The first alternative algorithm we considered used a massively-parallel array of line-segment finders (based e.g. on Xilinx gate arrays) that could process the resulting large number of combinations in $`<`$132 ns/crossing. The linefinders would look for hit triplets in three successive stations in each view, using linear interpolation since track curvature is negligible over the several-cm 3-station separation. This approach might be viewed as replacing detector slicing in $`\varphi `$ with slicing in $`z`$: a separate linefinder would be provided for each station, all running in parallel, and a subsequent step would link the track segments together into tracks, find the primary vertices, and identify any tracks that miss the primary.

The two-plane algorithm just described would find each track many times, since it looks in all stations in parallel. We have also explored alternative approaches that reduce the total processing power needed by finding each track only once near the point of optimal vertex resolution. Such an algorithm Erik  is the one we now propose to use BTeV-proposal . In this approach, pixel chips closest to the beam are classified as “inner” and those farthest from the beam (near the periphery of each pixel plane) “outer,” and only those track segments are sought that include inner pixels or outer pixels (see Fig. 3). A track segment including an inner pixel (inner segment) is required to extrapolate into the beam hole at the $`z`$ of the previous station, and one including an outer pixel (outer segment) is required to extrapolate beyond the vertex detector at the $`z`$ of the next station. To suppress ghost tracks, each inner segment must be confirmed by a corresponding outer segment. Each track is thus found once, near its points of entry into and exit from the vertex detector, and finding the remaining hits of the track, which do not provide essential information at this stage (although they can be used to improve the precision of the vertex fit), is postponed until Level 2. This algorithm can be implemented in an array of about 500 large FPGAs. The algorithm has been simulated and finds $`>`$90% of tracks while keeping the number of ghost tracks to an acceptable level.

Following the FPGA trackfinding stage just described, tracks are grouped into primary (interaction) vertices in an array of processor farms, comprising about 2500 DSPs in our current baseline design. Tracks with large impact parameters with respect to these vertices are evidence of heavy-quark decays. To reject light-quark background, events are required to contain $`n`$ tracks missing a primary vertex by at least $`m`$ standard deviations. (To avoid excessive false triggers from multi-interaction crossings, only impact parameters $`<`$2 mm are accepted.) A requirement such as $`n=2,m=6`$ rejects 99% of light-quark events while maintaining efficiency well above 50% for most $`B`$ decays of interest. The trigger system will have sufficient flexibility to allow a variety of such $`n,m`$ requirements to be used simultaneously, each with its own prescale factor.

The Level 2 and 3 triggers can be implemented in farms of commercial general-purpose processing nodes; we anticipate about 2500 such nodes in Levels 2 and 3. The goal of Levels 2 and 3 is to provide an additional factor of $``$100 in data rate while degrading efficiency for events of interest only minimally ($`<`$10%). To accomplish this, most remaining background events will be rejected, with much of the raw data in the remaining events summarized to reduce event size.

## III Trigger-DAQ Integration

The data rates into and out of the trigger levels are large (Table 1) and require close integration between the trigger and data-acquisition systems. These rates are carried by about 2500 1 GByte/s data links into and out of about 1 TByte of buffer memory. Detectors all read out directly into the Level 1 buffer, which passes the pixel data to the Level 1 trigger processors. For events failing Level 1, all data are flushed, while for those events that pass, additional tracking data are then provided from the Level 1 buffers to the Level 2/3 farms. As events move into the later stages of Level 3, data from additional detectors are provided to allow the decision to be further refined. Finally, events satisfying all trigger levels are compressed and written to an archival medium.

## IV Conclusions

To achieve unprecedented sensitivity to rare effects in charm and beauty decay, the BTeV collaboration has proposed to build the most sophisticated trigger system yet attempted in a high-energy-physics experiment. While mutiple reviews have revealed no show-stoppers, we realize that the task is challenging. Even downloading, control, monitoring, and calibration will require the development of new techniques CHEP .