text:
# SACY - Present Status 1

## 1. Introduction

The detection of X-ray sources by the ROSAT All-Sky Survey (RASS) associated with TTS outside star formation regions (Neuhäuser 1997) gave a tool to find new young associations. In fact, Torres et al. (2000), and Zuckerman & Webb (2000), using these sources, found evidences for two young associations near the South Celestial Pole, the Horologium (HorA) and the Tucana (TucA) Associations. To examine the possibility that these associations are the same and to search for other ones we undertook a Search for Associations Containing Young-stars (SACY) (de la Reza et al. 2001; Torres et al. 2001; Quast et al. 2001). In the SACY we selected and observed all bright RASS sources that could be associated with TYCHO-2 or HIPPARCOS stars with (B-V) $`>`$ 0.6, excluding very well known RS CVn, W UMa, giants, etc in SIMBAD.

We obtained high resolution spectra for the selected candidates with FEROS echelle spectrograph (Kaufer et al. 1999) (resolution of 50000; spectral coverage of 5000 Å) of the 1.52 m ESO telescope at La Silla or with the coudé spectrograph (resolution of 9000; spectral coverage of 450 Å, centered at 6500 Å) of the 1.60 m telescope of the Observatório do Pico dos Dias. For some stars we obtained radial velocities with CORALIE at the Swiss Euler Telescope at ESO (Queloz et al. 2000).

¿From the collected spectra we have obtained spectral classifications, radial velocities and equivalent widths of Li I lines. In particular, the Li I line is important since it can provide a first age estimate (Jeffries 1995) for late type stars allowing us to select possible Post-T Tauri stars. For a given star, if its Li I resonance line equivalent width is located near or above the Li I line delimited by the members of the Local Associations clusters (Neuhäuser 1997), it is flagged as young star.

## 2. Results

There are 9574 ROSAT bright sources in the Southern Hemisphere, 2071 of them with $`BV`$$`>`$ 0.6 in TYCHO-2. Until the La Serena meeting (2002, March) we observed 687 stars, and for 170 more stars we got relevant data from the literature. ¿From all these stars, 283 fulfill our definition for a young star (239 from our spectroscopic observations).

As most of these young stars have no HIPPARCOS parallaxes, we applied the following kinematical analysis to find possible associations: Each point in UVW velocity space is taken as a convergence point and we calculated for it the parallaxes of all stars such as to minimize the moduli of the space velocity vectors relative to this point (but, of course, preserving the parallaxes of HIPPARCOS stars). Then we calculated the density of stars in the velocity space around each point of a grid in UVW. Around some points there are density concentrations much larger than the background fluctuation, revealing possible associations. All the main concentrations are also constricted in space, but some of them cover large areas in the sky.

## 3. The Young Associations

In Tables 1 & 2 we present the properties of the eight young associations we found until now. For those that engulf previously known ones, we use bonafide members (not observed in the SACY) to help in their definition. The quantity of bonafide members used are indicated as the last number in Table 1. In Table 2 we indicate the distance ($`\rho _{max}`$) of the farthest member to the calculated center of the association, giving an idea of its size. It should be noted that the method gives no unique solution for some stars, but in almost all cases there is an obvious better membership. As there are important areas not covered until now by the SACY, and specially the Sco-Cen Association region, these identifications are preliminary.

The young association found up to now are the following:

a) GAYA \- The Great Austral Young Association (Torres et al. 2001) seems to split in two parts, although it is not clear if both are actually distinct associations or if this split was created by some bias or by our method. Their separation is mainly in W velocity and in distance. The previous HorA and TucA are within GAYA. Some of the proposed members of TucA are outside of the velocities definition (mainly the eastern ones). Until now both GAYAs seem deficient in binaries. There is a proposed SB2 member for GAYA2 but a radial velocity search on 28 members (more than a third of the members) has give no other SBs. The western proposed members are near the present survey limit.

b) TWA \- Until now we have found no new members. Of the known ones only five fulfill the SACY criteria. Our kinematical solutions put some stars farther away. TWA12 and TWA19 would be at about 100pc.

c) YSSA \- There is a group of young stars, spread from $`\rho `$ Oph to R Cra, with very similar properties, that we are now calling the Young Sco-Sgr Association. As we have no data yet for the Upper Sco subgroup we can not be sure if YSSA belongs or not to it. The eastern border of YSSA is at one of the present limits of SACY. The western border engulfs the stars mentioned in Quast et al. (2001) and Neuhäuser et al. (2000). Note the very well defined distance, near the assumed for the R CrA clouds.

d) $`ϵ`$ ChaA \- This association is defined by Mamajek, Lawson & Feigelson (2000). We propose new members, enlarging it. The well defined distance found by us indicates it is in front of the Cha complex. As in YSSA, the northern border being in the Sacy’s limits, we can not be sure if this group is not a part of Lower Cen-Crux subgroup.

e) $`\beta `$ PicA \- As described by Zuckerman et al.(2001) this association is very close to the Sun, containing members with distances between 10 and 50 pc. We proposed new members, someones at larger distances, among them V4046 Sgr, a notorious object, classified before as an isolated SB Classical T Tauri star.

f) OctA \- Is a very homogeneous small group of almost aligned stars (all young G stars) near the South Celestial Pole. As this region belongs to a completely surveyed area of the SACY, new members have too be found by other means.

g) AnA \- Another small group, very concentrated in velocity but spread out in space. Having only two stars with measured parallaxes, it should be viewed with caution.

### Acknowledgments.

This work was partially supported by a CNPq - Brazil grant to C. A. O. Torres (pr. 200356/02-0)

## References

de la Reza, J. R., da Silva, L., Jilinski, E., Torres, C. A. O., Quast, G. R. 2001, in ASP Conf. Ser. Vol. 244, Young stars near earth: progress and prospects, ed. R. Jayawardhana & T. P. Greene (San Francisco: ASP), 37

Jeffries, R. D. 1995, MNRAS, 273, 559

Kaufer, A., Stahl, S., Tubbesing, S., Norregaard, P., Avila, G., Francois, P., Pasquini, L., & Pizzella, A. 1999, Messenger, 95, 8

Mamajek, E. E., Lawson, W. A., & Feigelson, E. D. 2000, ApJ, 544, 356

Neuhäuser, R. 1997, Science, 276, 1363

Neuhäuser, R., Walter, F. M., Covino, E., Alcalá, J. M., Wolk, S. J., Frink, S., Guillout, P., Sterzik, M. F., Comerón, F. 2000, A&AS, 146, 323

Quast, G. R., Torres, C. A. O., de la Reza, J. R., da Silva, L., Drake, N. 2001, in ASP Conf. Ser. Vol. 244, Young stars near earth: progress and prospects, ed. R. Jayawardhana & T. P. Greene (San Francisco: ASP), 49

Queloz, D., Mayor, M., Naef, D., Santos, N., Udry, S., Burnet, M., Confino, B. 2000, in VLT Opening Symposium ¿From Extrasolar Planets to Brown Dwarfs, ESO Astrophys. Symp., Springer Verlag, Heidelberg, 548

Torres, C. A. O., da Silva, L., Quast, G., de la Reza, R., & Jilinski, E. 2000, AJ, 120, 1410

Torres, C. A. O., Quast, G. R., de la Reza, J. R., da Silva, L., Melo, C. H. F. 2001, in ASP Conf. Ser. Vol. 244, Young stars near earth: progress and prospects, ed. R. Jayawardhana & T. P. Greene (San Francisco: ASP), 43

Zuckerman, B., Sing, I., Bessell, M. S., & Webb, R. A 2001, ApJ, 562, L87

Zuckerman, B., & Webb, R. A. 2000, ApJ, 535, 959


format: text

text:
# Electronic Structure of the 𝑌𝐻₃ Phase from Angle-Resolved Photoemission Spectroscopy

## Abstract

Yttrium can be loaded with hydrogen up to high concentrations causing dramatic structural and electronic changes of the host lattice. We report on angle-resolved photoemission experiments of the $`Y`$ trihydride phase. Most importantly, we find the absence of metal $`d`$-bands at the Fermi level and a set of flat, $`H`$-induced bands located at much higher binding energy than predicted, indicating an increased electron affinity at $`H`$ sites.

Recently, switchable optical properties of metal hydrides at ambient pressures and temperatures have attracted strong interest. For trivalent $`Y`$, for instance, up to three hydrogen atoms can be absorbed. The dihydride is even a better metal than $`Y`$ itself but during the transition to the trihydride phase it turns from shiny metallic to transparent and insulating. With hydrogen a proton and an electron are introduced to the metal host. This results in doping the host material. Different models have been proposed to explain this spectacular metal-insulator transition. However, the behavior of $`H`$ in such hydrides is still under debate.

State-of-the-art *ab-initio* local density approximation (LDA) calculations do not reproduce the optical gap necessary to explain the transparent state in the trihydride phase unless additional symmetry lowering is considered by displacing $`H`$ atoms away from positions given by the $`HoD_3`$ structure of $`YH_3`$.

Other models, based on strong electron correlations, have been proposed to explain the metal insulator transition. Ng et al.  studied the effect of correlations on the bandwidth of hydrogen induced states. Hydrogen is present in the form of $`H^{}`$ as determined from electromigration experiments , where one electron is taken from the metal host. The two electrons on $`H^{}`$ are strongly correlated and the essence of the result of Ng et al.  is that the opening of the band gap is due to a correlation-induced band narrowing.

The model of Eder et al.  is based on the observation that the radius of the hydrogen is very sensitive to the occupation number. The two electrons on $`H^{}`$ are correlated but with drastically different radii around the proton. This results in a so-called breathing mode and a local singlet-like bound state with one electron on the proton and the other on the neighboring metal orbitals. Already at the meanfield level this introduces a significant correction to the potential at the $`H`$ site, effectively increasing the electron affinity, lowering the hydrogen band and opening the gap. There is an interesting and appealing connection between this model and the Zhang-Rice singlet  in high temperature superconductors (HTc’s). The two electrons form a singlet analogous to holes in the HTc’s.

On the other hand, very recent so-called GW-calculations  demonstrate the formation of a sufficient gap to explain the metal-insulator transition without need of strong electron correlations. Rather, these calculations indicate that the gap opening is described as in normal semiconductors. For semiconductors LDA does not produce the correct gap whereas the self-energy corrections included in the GW-calculations are able to account for this deficiency.

Indeed, detailed angle-resolved photoemission (ARPES) experiments are needed to favor one or the other model. However, practically all previous work on metal hydrides has been done on polycrystals and/or on samples that are capped with a protective Pd layer. In order to perform ARPES experiments, uncapped single crystalline material is needed. Furthermore, preparation has to take place *in situ* since $`Y`$ is extremely reactive.

Here we present, to our knowledge, the first ARPES data on the trihydride phase of uncapped single crystalline films. We find that the overall bandwidth of the Y trihydride phase agrees with LDA calculations. However, a set of flat bands is observed with significantly higher binding energy, a fact that argues in favor of ideas of the model proposed by Eder et al. .

Experiments were performed in a Vacuum Generators ESCALAB Mk II spectrometer with motorized sequential angle-scanning data acquisition.  ARPES measurements were performed with monochromatized He I$`\alpha `$ ($`h\nu =`$21.2 eV). The setup including a plasma discharge lamp with monochromator is described elsewhere . The energy resolution is 35 meV for He I$`\alpha `$ measurements. All measurements have been done with the sample kept at room temperature. The substrate used to grow on the single crystalline $`Y`$ hydride films was a W(110) single crystal. Details of the set-up used, the preparation and the careful characterization of the crystal and electronic structure and the calibration of the $`H`$-concentration has been described elsewhere . In brief, this set-up allows to prepare clean, single-crystalline rare-earth hydride films using hydrogen pressures up to 1.3 bar. It combines a high-pressure reaction cell with a custom made hydrogen purification system based on a Pd-24%Ag permeation tube and a sorption pump, i.e., a getter alloy (70%Zr-24.6%V-5.4%Fe). The hydrogen purification system removes residual contaminations from 6N $`H_2`$ gas efficiently. The overall design is such that the sample never gets in contact with non ultra-high-vacuum compartments. The results presented here are taken from 200 Å thick, well-ordered single crystalline films in the trihydride phase. The hydrogen composition is determined via photoelectron diffraction and X-ray photoelectron spectroscopy (not shown) to be $`YH_{2.9}`$.

Bandstructure calculations using the full-potential linearized augmented plane-wave method  within the generalized-gradient approximation  have been performed for comparison with the experiments. Calculations were performed for the $`HoD_3`$ structure with space group 165 ($`P\overline{3}c1`$). The lattice parameters used for the $`YH_3`$ calculation were $`a=b=6.34`$Å and $`c=6.6`$Å. A total of 485 k-points within the irreducible wedge of the Brillouin zone (BZ) were considered for the self consistency cycles and convergence was reached to within 0.1 mRy.

Figure 1 shows ARPES data measured along two high symmetry directions of the surface BZ. The situation in reciprocal space is illustrated in Fig. 2 where the surface BZ (left) and a cut through the extended bulk Brillouin zone containing the $`\mathrm{\Gamma },A,H,L,M,K`$ high symmetry points (right) is shown. $`YH_3`$ with the $`HoD_3`$ structure consists of a $`(\sqrt{3}\times \sqrt{3})R30^\mathrm{o}`$ reconstruction with respect to the hexagonal $`Y`$-lattice, i.e., the unit cell has $`a`$ and $`b`$ vectors which are $`\sqrt{3}`$ times longer and rotated by $`30^\mathrm{o}`$. Therefore the $`YH_3`$ BZ is smaller and rotated accordingly (Fig. 2).

Inspecting the experimental spectra (Fig.1 ) a similar behavior is observed for both directions. The dispersion appears relatively weak and the electronic states do not reach the Fermi level ($`E_F`$) or $`0`$eV binding energy, which is indicative for a gap. The $`d`$-states that are present in $`Y`$ (\[Kr\]$`4d^15s^2`$) and $`YH_2`$ (not shown) have disappeared. Five dispersing features are easily discernible and are labeled from $`1`$ to $`5`$. A broad almost dispersionless maximum between 5.5 and 8.3 eV (shaded area) persists for all angles respectively $`k`$-points. At first sight one is tempted to interpret the fact that the dispersion is weak and the spectral features are broad as caused by strong electron correlations and self-energy effects.

However, taking a closer look reveals a different point of view. In Fig. 3 the result of our bandstructure calculation of $`YH_3`$ in the $`HoD_3`$ structure is displayed. As mentioned above the approach using LDA does not show a gap and the system is metallic. The calculation presented here is in good agreement with the published results of Kelly et al.. We notice that there are many bands. This is a consequence of the large real-space unit cell of the $`HoD_3`$ structure containing many atoms. In fact, there are 24 atoms including 6 yttrium atoms and 18 hydrogen atoms in the unit cell. The important point to notice here is that there is a set of rather flat, hydrogen-induced bands (as deduced from a band character analysis) between 1.8 and 4.6 eV (dark shaded area) extending over the whole bulk BZ. Such a set of flat bands naturally gives rise to a high density of states in this energy region and leads to the interpretation of the broad, high intensity features in the experiment (shaded area in Fig. 1) in terms of these flat bands. The energy width of the region in the experiment agrees well with the energy interval of the flat bands in the bandstructure calculation. The only difference is that the experimental bands occur at a 3.7 eV greater binding energy.

In order to obtain a more detailed comparison we use a gray scale representation of both, the experimental data and the calculated band structure. Figure 4 displays the results along the $`\overline{\mathrm{\Gamma }}\overline{K}`$-direction. For the experimental data (left) the second derivative of the spectra of Fig. 1 have been plotted as a function of the wave vector parallel to the surface. The reason for plotting the second derivative is to flatten the spectra and to accentuate the dispersing features. The calculation (right) follows the free electron final state wave vectors drawn in Fig. 2 . Again we can identify the flat bands (marked by shading the energy scales) in both the experiment and the calculation where they occur with a 3.7 eV smaller binding energy. In addition we can identify the bands labeled $`1`$ to $`5`$ as already indicated in Fig. 1. Therefore comparing our experiments with the LDA calculation shows that the overall band width compares well. The differences are a rigid shift of approximately 1.5 eV and the position of the flat bands. However, a detailed comparison of different bands and their relative positions depends on the validity of the free-electron final-state approximation. Therefore we cannot, at present, give a robust statement about the exact position of the top of valence band (label 1, Fig. 4, left) and, as a consequence, about the size of the gap. This depends on the exact location in the component of the wave vector perpendicular to the surface, $`k_{}`$. On the other hand, the difference in the position of the flat bands is independent of this approximation since the flat bands extend over the whole BZ as seen from experiment (Fig. 1) and theory (Fig. 3).

It is interesting now to compare our results with theoretical models. At first we look at the GW calculations. In the study of Miyake et al.  the gap opens due to a shifting up of the conduction bands whereas in the case of van Gelderen et al.  the gap opens because of both, shifting up of the conduction bands and shifting down of valence bands. Otherwise the quasiparticle wavefunctions are practically identical to the LDA wave functions and there is basically a one to one correspondence between the bands of the LDA and GW bandstructures. The total band widths do not differ significantly as well. From our experiments we conclude that the overall bandwidth is well reproduced by LDA and GW-calculations. However, the high binding energy of the flat bands is not reproduced. Even the downwards shift of the valence bands ($``$1 eV) in the case of van Gelderen et al.  is not sufficient and does not account for the $`relative`$ shift of certain bands. Nevertheless, from the comparison in Fig. 4 (therefore assuming the free-electron final-state to be valid) we notice that a downward shift of the LDA results, as would be obtained for the GW-calculation produces quite a good agreement if we do not consider the flat bands.

Clearly, the good agreement of the overall bandwidth with the LDA calculations is in contrast to the model of Ng et al.  where strong electron correlations induce a band narrowing. However, we cannot exclude that individual bands are subject to significant narrowing. The reason is that from the experimental data, basically due to the large number of bands, it is not possible to identify every individual band. The width of the experimentally observed spectral features is not resolution limited and might be explained either by self-energy effects or also by a relatively high mobility of hydrogen atoms and defects due to the substoichiometry of the films. Furthermore, there is still debate on whether the $`HoD_3`$ structure really corresponds to the one for $`YH_3`$.

Finally we compare our results to the model by Eder et al.. The dependence of the occupation number of electrons on the hydrogen site is incorporated in their model Hamiltonian by having a novel hopping integral. The resulting singlet-like bound state is manifest via increasing the binding energy of the hydrogen states. This model and the one by Ng et al.  have in common that they treat hydrogen as impurities in the background of the metal atoms. However, whereas in the treatment of Ng et al. the result is a narrowing of bands, Eder et al. predict a shift of the potential at the hydrogen site $`retaining`$ the broad hydrogen band.

It is now very tempting to interpret our data, showing a shift of bands towards higher binding energy, in terms of such an additional potential on hydrogen sites. Observing Fig. 3 it appears that only the set of flat bands is strongly shifted down whereas the others, covering the full bandwidth, are consistent with the GW approach. This might be interpreted in the sense that the bands within the set of flat bands have a more localized character and therefore are stronger affected by the singlet-like bound state. However, a detailed comparison of how different bands are affected is not possible since a realistic bandstructure for the $`HoD_3`$ structure is not available within Eder’s model Hamiltonian.

In conclusion we have performed ARPES experiments on single crystalline films of Y in the trihydride phase. A comparison with LDA calculations shows a good agreement of the overall bandwidth. No significant band narrowing is observed although a narrowing of individual bands cannot be excluded. However a set of hydrogen induced, flat bands are identified in the experiment with a binding energy $``$ 3.7 eV higher than in the calculation. Furthermore, based on the free-electron final-state approximation, we infer a rigid shift of the valence band towards higher binding energy by more than 1 eV which is consistent with GW-calculations  and the opening of a gap. The downwards shift of flat bands supports the model proposed by Eder et al. predicting a shift of the potential at the hydrogen site and $`retaining`$ the broad hydrogen band.

Skillful technical assistance was provided by E. Mooser, O. Raetzo, R. Schmid, O. Zosso, Ch. Neururer and F. Bourqui. This project has been supported by the Fonds National Suisse de la Recherche Scientifique.


text:
# Extreme Self-Organization in Networks Constructed from Gene Expression Data

(Received 31 July 2002; published 12 December 2002)

## Abstract

We study networks constructed from gene expression data obtained from many types of cancers. The networks are constructed by connecting vertices that belong to each others’ list of $`K`$-nearest-neighbors, with $`K`$ being an a priori selected non-negative integer. We introduce an order parameter for characterizing the homogeneity of the networks. On minimizing the order parameter with respect to $`K`$, degree distribution of the networks shows power-law behavior in the tails with an exponent of unity. Analysis of the eigenvalue spectrum of the networks confirms the presence of the power-law and small-world behavior. We discuss the significance of these findings in the context of evolutionary biological processes.

preprint: Phys. Rev. Lett. 89, 268702 (2002); cond-mat/0207409

Recent technical advancements have led to widespread use of gene chips for quantizing and monitoring the expression level of thousands of genes in parallel Land99Gerh99 . Presently, gene expression profiling has become an important tool for diagnosis and classification of diseases. It is being used extensively for identifying genes responsible for specific conditions, e.g., various cancers Alon99 ; Nott01 ; Golu99 ; Pero99 ; Pero00 . This is done using specialized clustering techniques developed in recent years Eis98 ; Blat96 ; Braz00 . Gene expression data can also be used for constructing networks of coexpressed and coregulated genes. Since proteins are the end product of gene expression, various types of protein networks Jeo01Mas02  and gene networks are directly related. Consequently, networks of coregulated genes are expected to play key role in biological processes. In this letter we outline results that show the relevance of these networks in evolutionary biological processes.

The volume of gene expression data obtained from typical experiments is enormous and contains information on expression of all the genes (presently almost 10 000 or more) marked on the chip. In any given condition most of the genes are not important and do not express. As a result, before the expression data can be used for constructing networks, it requires extensive processing and filtering to eliminate uninformative genes. We skip these details here as they can be found with the source of the data Alon99 ; Nott01 ; Golu99 ; Pero99 ; Pero00  and elsewhere Eis98 ; Braz00 . Henceforth, we assume that expression data for the selected set of informative genes is available in the form of a matrix having $`N`$ rows and $`D`$ columns. The rows represent the genes and the columns represent the samples/tissues. Furthermore, the expression values of each of the genes in this matrix are normalized to have a mean of zero and variance of unity across the samples.

The normalized expression levels are treated as coordinates of $`N`$ genes, that form the vertices of the networks, in $`D`$ dimensional space of samples. The network construction algorithm requires specification of the maximum number of neighbors $`K`$, $`0K<N`$, that a vertex can have. For a given $`K`$, gene network is constructed using the following two step procedure. (i) For each vertex $`i`$, $`i=1`$, …, $`N`$, make a list $`L_i`$ of its first $`K`$ nearest neighbors. (ii) Connect all vertices $`i`$ and $`j`$ through an edge if $`iL_j`$ and $`jL_i`$, otherwise the vertices are not connected. This algorithm is derived from the $`K`$-nearest-neighbor parameter estimation method Fuk90 . With some heuristic modifications it has been used earlier in clustering analysis of various types of data Blat96 . We used Euclidean norm as the distance measure for making the list of $`K`$ nearest neighbors. Other distance measures can also be used. The results presented herein remain unaltered as long as the distance measure preserves the ordering of points obtained from the Euclidean measure.

For a given data set, the topological structure of networks generated by this algorithm depends strongly on the parameter $`K`$. For $`K=0`$, the network consists of $`N`$ isolated vertices, and for $`K=N1`$ each vertex is connected to all the other vertices. For most of the values of $`K`$, these networks have more than one connected component. As $`K`$ increases, the connectivity of each vertex grows depending on its local environment. Vertices that lie close to each other tend to get mutually connected in preference to those lying farther away. Thus, all connected components in these networks have a small-world structure. Furthermore, for $`K3`$ the networks have a giant connected component.

We analyzed networks constructed using several published gene expression data sets. To ensure that our results are not affected by possible bias of technology used for manufacturing the gene chips, we used expression data sets obtained using both oligonucleotide arrays Alon99 ; Nott01 ; Golu99  as well as cDNA microarrays Pero99 ; Pero00 . For each network we calculated $`P(z)`$ the probability of finding a vertex of degree $`z`$. $`P(z)`$ is normalized by $`N`$ so that $`_zP(z)=1`$. Since these networks are small, $`P(z)`$ is very noisy. As a result, we calculated the cumulative probability distribution

$$F(z)=\underset{i=z}{\overset{z_{\text{max}}}{}}P(i)$$
(1)
where $`z_{\text{max}}`$ is the maximum degree in the network.

We also define a quantity $`c=z+1`$. This quantity gives the size of the smallest cluster around a vertex with connectivity $`z`$ that includes the vertex and its neighbors. It can also be considered as the size of a “droplet” that is formed by the vertex and its neighbors. Since, the smallest and the largest values of $`z`$ are 0 and $`z_{\text{max}}`$, the corresponding values for $`c`$ become 1 and $`c_{\text{max}}=1+z_{\text{max}}`$, where $`c_{\text{max}}`$ is the largest droplet size. The probability density $`\stackrel{~}{P}(c)`$ and cumulative distribution $`\stackrel{~}{F}(c)`$ corresponding to $`c`$ are defined similarly to those for $`z`$. The linear relationship of $`c`$ and $`z`$ implies that $`\stackrel{~}{P}(c)=P(z)`$ and $`\stackrel{~}{F}(c)=F(z)`$. Outside the range of $`z`$ and $`c`$, corresponding probability density functions are defined to be zero. Thus, $`\stackrel{~}{F}(c)|_{c1}=F(z)|_{z0}=1`$ and $`\stackrel{~}{F}(c)|_{c>c_{\text{max}}}=F(z)|_{z>z_{\text{max}}}=0`$.

The homogeneity of the networks can be characterized by a single order parameter. A suitable candidate for this is the area $`\mathrm{\Lambda }`$ enclosed by the $`\stackrel{~}{F}(c)`$ versus $`c^{}=c/c_{\text{max}}`$ curve between $`c^{}=0`$ and $`1`$. It takes values in the range $`0`$ to $`1`$ depending on the homogeneity of the network. Since $`c_{\text{max}}`$ is finite and the values of $`c`$ are evenly spaced, $`\mathrm{\Lambda }(K)`$ is easily calculated using Trapezoidal rule and equals

$$c_{\text{max}}\mathrm{\Lambda }(K)=\frac{1}{2}\left[1P(z_{\text{max}})\right]+\underset{i=0}{\overset{z_{\text{max}}}{}}(i+1)P(i).$$
(2)
The terms with factor of $`1/2`$ represent the area of strips at the boundary. Their contribution vanishes as the number of strips increases. It is zero if $`P(z)`$ is a delta function.

From Eq. (2) the value of $`c_{\text{max}}\mathrm{\Lambda }(K)`$ is easily identified as the mean size $`\overline{c}=1+\overline{z}`$ of the droplets, where $`\overline{z}`$ is the average degree of vertices in the network. Thus, $`\mathrm{\Lambda }(K)`$ is the average droplet size normalized with the size of the largest droplet in the network. It measures separation between the mean and maximum droplet sizes in the network and functions as an indicator of the the overall behavior of $`F(z)`$ and $`P(z)`$. Very small values of $`\mathrm{\Lambda }(K)`$ imply that $`F(z)`$ descends sharply from unity to almost zero at a small $`zz_{\text{max}}`$ and becomes nearly flat with plateau stretching until $`z_{\text{max}}`$. The corresponding $`P(z)`$ has long tails with weight centered at small $`z`$. $`\mathrm{\Lambda }(K)`$ close to unity implies that $`F(z)`$ stays nearly flat at unity for most of $`zz_{\text{max}}`$ and descends sharply to zero at some $`zz_{\text{max}}`$. In this case $`P(z)`$ is sharply peaked similar to a delta function near $`z_{\text{max}}`$. Intermediate values of $`\mathrm{\Lambda }(K)`$ imply a decaying $`F(z)`$ corresponding to various forms of $`P(z)`$ including those that rise slowly to a peak and then decay slowly via sharply truncated power-law tails.

Figure 1 (Alon99) shows the behavior of $`\mathrm{\Lambda }(K)`$ for networks constructed using colon cancer data of Alon et al. Alon99 . The figure shows, that as $`K`$ is increased from 1 to $`N1`$, $`\mathrm{\Lambda }(K)`$ initially decreases and attains a minimum at some value of $`K=K_1`$ (here $`K_116`$). The minimum is nearly flat and persists until $`K=K_2`$ (here $`K_224`$). As $`K`$ is increased beyond $`K_2`$, $`\mathrm{\Lambda }(K)`$ continues to increase, reaching its maximum value of 1 at $`K=N1`$. Fig. 1 (Pero99-26) shows similar behavior of $`\mathrm{\Lambda }(K)`$ in networks constructed using breast cancer data of Perou et al. Pero99  with a second minimum at $`K3N/4`$. The second minimum is usually a finite size effect. It can also occur if there is high inhomogeneity on length scales large compared to the nearest neighbors scale. It is very shallow compared to that between $`K_1`$ and $`K_2`$ because the probability density function corresponding to the size of “droplets” on larger length scales are sharper compared to $`P(z)`$ and also have smaller tails. The higher order minima, however, are not relevant. Behavior similar to that seen in Fig. 1 was observed in networks from several other gene expression data sets also Nott01 ; Golu99 ; Pero00 .

The behavior of $`\mathrm{\Lambda }(K)`$ divides the networks in three classes. (i) The networks corresponding to $`1K<K_1`$ have few connections between vertices but have high homogeneity. (ii) The networks for $`K_1KK_2`$ are somewhat better connected and highly inhomogeneous. (iii) In the networks for $`K_2<KN1`$, the vertices have many connections and high homogeneity that increases with $`K`$.

Figure 2 shows the variation of the observed cumulative probability distribution $`F(z)`$ with the normalized degree $`(z+1)/(z_{\text{max}}+1)`$ in a wide range of values of $`K`$ for networks constructed from many gene expression data sets Alon99 ; Nott01 ; Golu99 ; Pero99 ; Pero00 . It is clear from the figure that all the curves in the range $`K_1KK_2`$ (solid circles) show a very good collapse, and thus good scaling, for all the data sets. This range of $`K`$ houses the minimum of the order parameter and the corresponding networks are highly inhomogeneous with loosely connected vertices. The solid straight line passing through the tails of these curves is a least-square fit of the form

$$F(z)=ab\mathrm{ln}\left(\frac{z+1}{z_{\text{max}}+1}\right)$$
(3)
where $`a`$ and $`b`$ are the fit parameters. For the graphs shown in the figure, these parameters vary in the range $`0.03a0.08`$ and $`0.47b0.71`$ over all the data sets.

The extremely good fit of Eq. (3) in the tails of $`F(z)`$ seen in Fig. 2 implies that the corresponding probability density functions have a scale-free behavior of the form $`P(z)b(z+1)^1`$ in the tails. The power law is seen in the range $`0.35(z+1)/(z_{\text{max}}+1)1`$, i.e., in nearly 60% to 65% of the range of variation of vertex degree. As this range is small for observing heavy tailed distributions, we analyze the eigenvalue spectrum of the adjacency matrix of networks Far01Goh01  to have another evidence of the scale-free character of the networks.

The spectral density $`\rho (\lambda )`$ of the eigenvalue spectrum of the adjacency matrix of networks

$$\rho (\lambda )=\frac{1}{N}\underset{j=1}{\overset{N}{}}\delta (\lambda \lambda _j),$$
(4)
where $`\lambda _j`$ is the $`j`$th eigenvalue of the adjacency matrix, is a good indicator of the overall behavior of their degree distribution $`P(z)`$ and topological structure. For random graphs having a giant connected component $`\rho (\lambda )`$ is known to converge to a semicircle following the Wigner’s law Far01Goh01 . Deviations from Wigner’s law are seen for other cases. For the small-world networks $`\rho (\lambda )`$ shows a complex highly skewed structure with several blurred peaks Far01Goh01 . For scale free networks, the spectral density has a triangular shape with the central part lying above the semicircle Far01Goh01 .

Figure 3 shows the spectral density of networks constructed from colon cancer data of Alon et al. Alon99  for different values of $`K`$ corresponding to the three zones of behavior of $`\mathrm{\Lambda }(K)`$ seen previously. The figure shows that for small $`K`$ the spectral density has an irregular shape with several blurred peaks and its bulk is confined below the semicircle. This is a characteristic of small-world networks Far01Goh01 . As $`K`$ is increased, the bulk portion starts assuming a triangular shape which is a little skewed for small $`K`$. The top portion of the triangle starts protruding above the semicircle as $`K`$ crosses $`K_1`$. This is clear from the intermediate and high values of $`K`$ in the figure. The triangular shape persists till $`K`$ becomes almost equal to $`N1`$. At $`K=N1`$, the spectral density develops a delta function peak corresponding to $`N1`$ repeated eigenvalues at $`\lambda =1`$ and the largest eigenvalue is $`\lambda _1=N1`$. This behavior of the spectral density was seen in all the other data sets also Nott01 ; Golu99 ; Pero99 ; Pero00 . It confirms that these gene networks have small-world character and become scale-free for $`KK_1`$.

Presence of scale-free behavior indicates high degree of self-organization in the system and is known to be a characteristic of natural systems Bak87 . It has been observed in several natural and artificial networks, e.g., power grid Wat98 , the Internet Fal99Coh00 , the World Wide Web Alb99Hub99 , actor network Bar99 , web of human sexual contacts Lil01 , citation and collaboration networks Red98New01 , conceptual network of language Mot02 , metabolic network Jeo00 , food web Sol00 , and protein networks Jeo01Mas02 . The exponents of the power laws observed earlier, however, were always more than unity. The gene networks studied here are first examples of biological networks showing scale-free behavior with exponent of unity.

The scale-free character of coexpressed gene networks means that these networks are extremely inhomogeneous and contain few genes that are very highly connected and a large number of genes with low connectivity. This implies that these networks contain large groups of coexpressed genes. As a result, the present study conclusively shows, using direct experimental data, that various types of cancers are a consequence of malfunction of only a few genes that either regulate the expression of a large number of other genes or form the hubs of various crowded gene regulatory pathways functioning in the organism. It is known that disturbing such genes could be fatal for the organism Has01 , which also turns out to be the mechanism of origin of cancers. Identification of such genes and understanding their functionality under various conditions through which an organism can pass in its lifetime is directly relevant in the design of highly targeted drugs, among other possibilities.

The simultaneous presence of small-world and scale-free characters in these networks seems to be a perfect fit in the evolutionary scheme of biological systems. The high robustness displayed by biological systems is a consequence of the scale-free character of the associated networks. On the other hand, the fast reaction and rapid adaptability shown by biological systems can come only if the associated networks have a small-world character. This fits the structure of biological signaling system well because the chemical signaling employed at most places in biological systems, by its very nature, is very slow compared to, e.g., electrical signaling in neurons. Thus, for achieving fast message transmission, the associated networks must evolve to have small-world character.

I thank Eytan Domany for introducing me to gene expression data analysis and helpful discussions, Itai Kela for help in filtering breast cancer data, and Deepak Dhar for critical review of the manuscript.


text:
# CLPGUI: a generic graphical user interface for constraint logic programming over finite domains

## 1 Introduction

Several tools for visualizing the execution of constraint programs have been developed in the last few years. These tools have been found very useful for debugging and improving constraint programs, and for teaching constraint programming. One can distinguish:

* post-mortem visualization tools, these tools are used after execution of the program, the program is annotated with specifications of the information to trace. This approach is implemented for example in the CHIP or CIAO systems, it allows using a wide variety of viewers, including both application oriented tools , and generic tools, for visualizing the search tree , finite domain variables , or constraint propagation .
* dynamic visualization tools, these tools are connected to the constraint programming interpreter and realize a visualization on-line, possibly with animations . This approach is implemented in Grace tool  for finite domains visualization, and in OPL studio  for search tree and constraint propagation visualization.
* dynamic visualization and control tools which allow interaction with a CLP process through different visualizations. One example is the Oz-Explorer system  where it is possible to jump to any previously encountered state by simply clicking on a node of the search tree, and restart computation from that state. User-guided search is implemented in Oz-Explorer using the first-class computation spaces of Oz. Recomputation is used to trade space for time in Oz-Explorer, and similarly in OPL studio , the state restoration mechanisms in tree search are described in .

In this paper we propose to push forwards these ideas towards an open architecture for connecting a CLP process with dynamic visualization and control tools. To our knowledge, the visualization of search trees in three dimensions has not been much investigated. We argue that the 3D representation of search trees proposed in this paper provides the most appropriate visualization of large search trees.

Our ambition is not to realize an ad hoc tool limited to a particular constraint programming system, namely GNU-Prolog , but a generic tool which can be ported to other constraint programming systems. The main reason for this is that a wide variety of viewers can be useful for debugging constraint programs and it is possible in this way to share developments. Our approach relies on the generic trace format which is defined in the OaDymPPac consortium  for post-mortem analysis. We propose to extend this format with a similar generic format for control, and to use these formats for connecting on-line the CLP process to the GUI process. Our implementation of CLPGUI  in GNU-Prolog and Java is depicted in Figure 1.

In Section 2 of this paper, we describe the client-server architecture of CLPGUI and the user console which facilitates the establishment of a connection with a CLP server, and controls the execution of the CLP program. Section 3 describes annotations that can be added to a CLP program, in order to give an external name to finite domain variables, to create buttons for posting constraints or goals from the CLPGUI console, and to define the level of granularity of the search tree which will be visualized. Section 4 presents a dynamic 3D viewer for visualizing the evolution of finite domains variables over time. Section 5 describes the representation of the explored search space by partial CSLD derivation trees, presents different visualizations with 2D and 3D viewers, and discusses their use in CLPGUI. Then in Section 6 we discuss the implementation of the interactive execution model and of the annotations, using the global variables of GNU-Prolog. Finally we provide some evaluation results, and conclude on the current limitations of the system and on some perspectives for future work.

## 2 Client-Server Architecture

The graphical user interface of CLPGUI is a Java application connected by sockets as a client to a server which executes CLP goals. Both processes can run on different machines and communicate over the network. This has been experienced with CLPGUI for visualizing the execution of CLP programs on a Workbench of Virtual Reality. The choice of the Java language for implementing the GUI is motivated by several reasons:

* its object-orientation, all 3D viewers presented in the following sections inherit from a single class for moving and projecting 3D figures;
* the encapsulation of events handling, that is preponderant in dynamic visualization;
* the threaded execution, which is mandatory for implementing communication with the CLP process;
* its wide availability.

For efficiency reasons, we did not use the Java-3D library for the viewers presented in this paper, as they can directly benefit from ad hoc optimizations that speed-up their incremental display. Nevertheless the architecture can support the use of the powerful Java-3D library for developing complex application-oriented viewers.

During initialization, the CLP server starts an interpreter of the command lines received on the socket. The command lines may contain any CLP goal. The GUI Java client opens a graphical console such as the one in Figure 2. That console is used for establishing a connection to the CLP server, and for posting constraints or executing CLP goals.

The CLP program may contain annotations for creating buttons for some constraints or for some Prolog goals to execute in an interactive manner. These buttons for posting constraints or Prolog goals then appear at the bottom in the CLPGUI console, see Figure 2. A click on the button posts the constraint or executes the goal associated to the button. Other arbitrary goals can be executed by entering them in a text field. In addition, one button called “backtrack” continues the execution of the current goal up to the next success, or, if there are no more success, returns to the state of the previous interaction. Another button called “backtrack interaction” forces backtracking to the state of the previous interaction. The menu bar of this console contains menus to select and activate the viewers of the search tree or of the finite domain variables.

## 3 Annotations

The CLP program may contain annotations for giving an external name to CLP(FD) variables, for creating buttons for posting constraints or goals from the CLPGUI console, and for marking the goals to visualize in the search tree. The following predicates are part of the annotation library:

* varnames(\[V1,...,VN\],\[name1,...,nameN\]) and varnames(\[V1,...,Vn\]) give an external name to a list of CLP variables. These external names are used in the graphical user interface and for the communication by sockets.
* button(goal) creates a button in the GUI console for executing a goal or for posting a constraint.
* bagof\_buttons(goal, call) creates a bag of buttons for each successful instance of the second argument.
* trace\_search(goal) executes the goal and traces the execution of that goal, by creating nodes in the search tree. The goals and constraints posted from the graphical console are always traced.
* show\_domains updates the visualization of the current state of FD variables.

Annotations provide a simple mechanism for encapsulating communications towards the GUI process . The advantages of this abstraction are:

* the flexibility for defining different levels of granularity concerning the information to visualize.
* the easiness for making interactive already existing programs,
* the portability of the GUI to other constraint programming systems, as all communications with the GUI are encapsulated in the implementation of annotation predicates.

The limitations of annotations are well-known in standard programming environments: they may be difficult to maintain in large programs. In that case, one solution is to automatically generate annotations with a graphical editor of the program source, where spy points and trace options can be specified. Nevertheless, one peculiarity of constraint logic programming is the conciseness of programs. CLP(FD) programs for solving combinatorial optimization problems on real-size data may compute with a huge amount of constraints and variables, but the program source for handling constraints and defining complex search strategy usually remains relatively concise. Therefore in this context, the proposed annotations appear as a satisfactory solution.

In many CLP systems however, the heuristic labeling procedures are built-in, and cannot be precisely traced with a simple annotation. In these cases, the annotations have to rely on the tracing facilities of the CLP system in order to extract, and communicate to the GUI, the relevant information. A simpler solution is to program the labeling heuristics in the host language, for making available in the host language the information coming from the constraint solvers that is relevant to the search heuristics. In that case, the effect of the search strategy can be visualized at different levels of granularity. In its simplest form, a predicate for tracing a labeling procedure can be defined with a trace\_search annotation as follows:

```
trace_labeling([]).
trace_labeling([X|L]):- trace_search(fd_labeling(X)), trace_labeling(L).

```

It is worth noting that even if the search strategy is implemented with a meta-interpreter, or constraint posting is combined with labeling, the relevant goals of the execution can still be traced with annotations.

The same difficulty arises for tracing internal constraint propagation steps. This is not possible without access to the wakening events of the constraint solver. The annotations for tracing constraint propagation steps have thus to rely on the tracing facilities of the solver in order to extract and communicate constraint wakening events.

###### Example 1

The following annotated GNU-Prolog program solves the well known SEND+MORE=MONEY puzzle in an interactive manner, by creating buttons for posting the constraints and for trying two labeling goals in this example:

```
sendmore(L):-
         L=[S,E,N,D,M,O,R,Y],
         varnames(L,[’S’,’E’,’N’,’D’,’M’,’O’,’R’,’Y’]),
         fd_domain(L,0,9),
         button(fd_domain([S,M],1,9)),
         button(1000*S+100*E+10*N+D + 1000*M+100*O+10*R+E
                #= 10000*M+1000*O+100*N+10*E+Y),
         button(fd_all_different(L)),
         button(trace_labeling(L)),
         reverse(L,L2),
         button(trace_labeling(L2)).

```

This program generates the console in Figure 2. The evolution of the finite domain variables over time, after the posting of constraints and of the first labeling goal, is depicted in Figure 4. The visualization of the search tree for obtaining all solutions under the first labeling goal, and then under the second labeling goal executed after a backtracking command, is depicted in Figure 3. Under the first ordering, the labeling is deterministic. Under the second ordering, few backtracking steps occur on variable Y when searching for other solutions. Note that other labeling heuristics can be tried directly from the console. On such pedagogical examples, the advantage of immediately visualizing the effect of posting a constraint or trying a labeling, is clear for teaching purposes.

###### Example 2

The following program solves the N queens problem by creating buttons for posting the constraints (safe predicate) and labeling goals for each variable (fd\_labeling predicate) and for all variables (trace\_labeling predicate).

```
queens(N,L):-
        length(L,N),
        varnames(L),
        fd_domain(L,1,N),
        button(safe(L)),
        bagof_buttons(fd_labeling(X),member(X,L)),
        button(trace_labeling(L)).

```

Three visualizations of the search tree for the 8 queens problem are depicted in Figures 5, 6 and 7.

## 4 3D Views of Finite Domain Variables

The evolution of finite domain variables is visualized in a three dimensional graph variable-domain-time, as proposed in the VIFID/TRIFID tool . Here the visualization is dynamic, the Java process reads the stream of finite domains information and paints the figure in an incremental manner. Domains are depicted by their size on the vertical axis, see Figure 4. According to options, that can be set in the CLP program or in the GUI, only the size, the interval or the complete domain of variables is visualized. But in any case only the sizes of the domains are memorized, therefore the extra information is lost when the figure is repaint. The time axis traces the interactions (i.e. the posting of constraints in the example), and the execution of traced goals (i.e. the labeling in the example). This view shows that the posting of constraints instantiate variables $`S,M,O`$ and that the first labeling step on variable E in fact instantiates all variables by constraint propagation. An option determines whether backtracked states are traced or erased.

The figure can be moved, zoomed and rotated. For efficiency reasons, the rotations are limited to a quadrant of a sphere which is not a real limitation for the user. In this way the visible faces are efficiently determined and the figure can be drawn incrementally.

Extra information on variables and executed goals can be obtained by moving the mouse on the position of a variable or on a time position.

The 3D dynamic view of finite domain variables evolution is very useful for teaching constraint programming. The effect of constraints is immediately seen and many strategies can be tried step by step. On large set of variables, the 3D view of domains can still be useful to get a view of the pruning power of different constraint modelings, and of the efficiency of different search heuristics, by comparing the general shape of domain reductions.

## 5 2D and 3D Views of the Search Tree

### 5.1 Partial CSLD derivation trees

The search tree considered in CLPGUI is a labeled tree defined as follows:

* a node is introduced for each call to a traced goal (called a call node), and for each success to a traced goal (called a success node),
* the label of a call node is the called goal,
* the label of a success node is the list of named variables with their value,
* the arcs correspond to the operational CLP transitions.

This tree is a subtree of the CSLD derivation tree . It is thus a quite natural representation of the search tree for describing CLP program execution. A branch represents a conjunction, and the different successors of a node represents a disjunction. A success node may have several successors if there is an untraced non-deterministic goal which is executed after the success, and before the next call to a traced goal. This is the main reason why success nodes are introduced in partial CSLD trees. In this way, the non-determinism due to untraced goals cannot be confused with the non-determinism of traced goals.

One disadvantage of CSLD trees is that in the case of deterministic programs they are threadlike and thus space consuming in their standard representation. AND-OR trees provide a more compact representation, as the threadlike parts of the CSLD tree are compacted in the successors of a single AND-node. For this reason, in the context of logic programs where most predicates are deterministic, AND-OR trees, and their variant AORTA diagrams which indicate the status of resolution of the goals, have been preferred . Nevertheless in the context of constraint logic programming over finite domains, the situation is quite different. The search tree to visualize is usually focused on the labeling predicates, or more generally on the branching procedure, which is highly non-deterministic (at least during debugging). The representation of the deterministic part of the search tree with threadlike structures provides an immediate visualization of the pruning power of constraints.

A naive solution for tracing constraint propagation steps in this approach is to add deterministic nodes for tracing constraint wakening events. For space limitation reasons, it is preferable however to aggregate constraint propagation information to the nodes of the search tree. This is proposed in the “Christmas trees” of OPL studio .

For search engines not based on backtracking, it is worth noting that a partial CSLD derivation tree can still provide a valid representation of the explored search space, as long as the explored states can be defined by their relation to some ancestor states. A formalization of an interactive constraint solver by transformations of CSLD derivation trees was done in .

### 5.2 Visualization

Once the search tree is formally defined, it can still be visualized in many ways, and in some cases it can be interesting to use several visualizations at the same time. We have currently implemented several two-dimensional and three-dimensional viewers, but many more representations could be imagined and fruitfully used.

In all the following representations, the labels of the nodes are visualized when the mouse is moved on them, and an option makes all nodes visible. The successes are materialized by a red cross. Each view can be moved, zoomed and rotated.

Figure 3 uses a standard 2D representation of the search tree in a fixed width. Figure 5 uses a dynamic 2D representation of the tree with a fixed spacing between leaves. This representation of the tree can be drawn incrementally and is thus appropriate for the dynamic visualization of large trees.

To our knowledge, the 3D visualization of search trees has not been much investigated. Figure 6 shows a somewhat original 3D representation of the search tree with alternating planes of successors. One advantage of this 3D representation is that it is relatively compact, it helps visualizing rather large trees by playing with rotations, see Figure 8 for another example. Our experience is that the 3D view is the most appropriate view to apprehend the shape of large search trees.

It is interesting to note that one obtains a treemap representation of the tree by rotation of the 3D alternate tree up to its vertical projection, as done in Figure 7. Treemap representations (with colors for aggregating information) are known to be particularly efficient to represent very large data  and to visualize complex phenomenons such as correlations, patterns or symmetries.

The interaction allowed in these views to restore a state is currently limited to user-guided backtracking and recomputation. The automatic recomputation of any state of the tree as described in  is currently not implemented.

### 5.3 Branch and bound optimization

The branch and bound procedure is widely used in constraint programming to solve optimization problems. Branch and bound optimization develops search trees in two parts. The first part corresponds to the enumeration of solutions with decreasing costs (for minimization problems). The second part exhausts the search space to show that there does not exist a better solution than the last solution found. The second part of the search tree constitutes the proof of optimality.

Figure 8 shows the search tree for the bridge problem , a medium size job-shop scheduling problem. The first descent corresponds to the search of the first solution of cost 108. It contains 78 call nodes. The second descent (after some hesitation) corresponds to the search of the optimal solution of cost 104. It contains 99 call nodes. The bottom part of the tree contains 3728 call nodes and corresponds to the proof of optimality. The 3D view is the most appropriate view of this search tree. It can be moved and rotated without difficulty.

## 6 Implementation

### 6.1 Interactive Execution Model for CLP

The interactive execution model of the CLP process derives from a more general model for adding and removing constraints and goals described in . In CLPGUI, constraints and goals can only be added to the current goal, the removing of constraints or goals occurs by backtracking. It is therefore possible on a success of the current goal:

* to add constraints or any goals to the current goal and continue resolution,
* to backtrack to the next success (command “backtrack” of Section 2),
* or to backtrack to the last interaction (command “backtrack interaction”).

It is worth noting that such a top level is in fact very appropriate for standard Prolog systems, where the capability of adding goals on a success of the current goal, and continue resolution, is usually missing.

Our current implementation uses the global variables of GNU-Prolog  to memorize global information, such as input and output sockets, variable names, and information used for backtracking. Global variables make it possible to avoid adding parameters to many predicates and lead to a simple implementation of annotations.

### 6.2 Communication messages

In this section we describe the communication messages which are transmitted between the CLP process and the GUI process. The CLP process produces the trace information specified by the annotations in the CLP program, or asked from the GUI. The messages emitted from the CLP to the GUI are the following:

* \<variables ...\> sends the list of FD variable names
* \<button G\> asks the GUI to create a button for posting the constraint or goal G
* \<undo button G\> indicates backtracking on the creation of a button for G
* \<node G\> traces a call to goal G
* \<undo node G\> traces backtracking on the call to G
* \<child G\> traces a success to G
* \<undo child G\> traces backtracking on the success to G
* \<undo goal G\> indicates backtracking on the call to goal G
* \<domainSizes ...\> sends the domain sizes of FD variables
* \<domainIntervals ...\> sends the current intervals of FD variables
* \<domainValues ...\> sends the current finite domains of FD variables
* \<undo domainValues\>, \<undo domainIntervals\>, \<undo domainSizes\> warns the GUI that the finite domain variables are updated by backtracking.
* \<success\> indicates that the current derivation is a success
* \<clear\> indicates return to top level.

In the other direction, the messages emitted from the GUI to the CLP process are the following:

* \<showSize\>, \<showInterval\>, \<showValues\> sets the information on finite domains that need be sent
* \<execute G\> asks to post constraint G or execute goal G
* \<backtrack\> asks backtracking to the next success
* \<backtrackInteraction\> forces backtracking to the last interaction
* \<clear\> asks to abort the current execution.

The portability of CLPGUI to a new constraint programming system is determined by the ability of the constraint programming system to produce and interpret these communication messages. The messages of the first list are produced by the predicates of annotation library described in Section 3. The messages in the second list are interpreted by the interactive execution model described in the previous section.

## 7 Evaluation

Our experience of using CLPGUI for teaching constraint programming has been very positive. The dynamic visualization of CLP programs really speeds-up the process of learning the basic concepts of domain filtering, constraint propagation and search trees. CLPGUI has also been fruitfully used to visualize the search tree of CLP(R) programs.

On real-size data, CLPGUI has shown satisfactory performance figures. We report in this section the timings on a Pentium III 600 MHz processor under Linux. GNU-Prolog solves the bridge problem mentioned in Section 5.3 in 100 ms, including the proof of optimality. The solving together with the visualization of the search tree with 3905 call nodes takes 470 ms with CLPGUI. This overhead is due to the communication of messages by sockets. The overhead was reduced from 2600 ms to 470 ms by optimizing socket calls and by using simple data compression techniques for communication. Moreover, the drawing of the tree is immediate and the figure can be moved and rotated without difficulty.

## 8 Conclusion and future work

We have described an open architecture for visualizing and controlling the execution of constraint logic programs. Communication by sockets between the CLP process and the GUI process has proved efficient enough for dynamic visualization and interactions. An important reduction of the overhead was obtained by optimizing socket calls and by using simple data compression techniques for communication.

CLPGUI supports the use of different viewers. Our experience has shown that a somewhat original 3D visualization of the search tree proposed in the paper, is often the preferred view to apprehend the shape of large search trees, as it is very compact. More work is needed however to parametrize the different viewers and invent novel visualizations of complex data. In this respect the flexibility of the architecture makes it possible to connect CLPGUI to external generic viewers, or to use powerful libraries like Java 3D to develop application-oriented viewers.

The most obvious limitation of our current implementation of CLPGUI is the absence of connection to a system for tracing constraint propagation, simply because such a tracer does not exist yet for GNU-Prolog. Nevertheless a generic trace format for finite domain constraint solvers has been defined in the OaDymPpac consortium , and we plan to use this format in future versions of CLPGUI.

Finally, CLPGUI is not a visual programming tool as far as the capabilities of defining goals from the GUI are extremely rudimentary. Nevertheless the proposed architecture can support this kind of extension by adding the capability to define constraints and goals graphically, that is certainly worth investigating.

#### Acknowledgements.

I would like to thank all members of the OaDymPpac RNTL project of the French Ministry of Research, and especially Abderrahmane Aggoun, Thomas Baudel, Pierre Deransart, Jean-Daniel Fekete, Ludovic Langevine and Mohammad Gonhiem for interesting discussions on this topic. I am grateful also to Anupam Agarwal for optimizing communication by sockets, and to Jean-Michel Leconte for his preliminary work along these lines on the workbench of virtual reality at INRIA.


text:
# 1 Introduction

## 1 Introduction

This paper presents the results of two types of search for the production of a di-photon system recoiling from another massive scalar or vector object. The searches are sensitive to the processes $`\mathrm{e}^+\mathrm{e}^{}\mathrm{XY}`$, with $`\mathrm{X}\gamma \gamma `$ and $`\mathrm{Y}\mathrm{f}\overline{\mathrm{f}}`$, where $`\mathrm{f}\overline{\mathrm{f}}`$ is a hadronic system (jets), a pair of charged leptons, or neutrinos resulting in missing energy. In the general search mode $`\mathrm{X}`$ must be a scalar, $`\mathrm{Y}`$ can be any scalar or vector particle of any mass, and both particles must be short-lived so that they decay close to the interaction point. The other search mode is referred to as the $`\mathrm{h}^0`$$`\mathrm{Z}^0`$ search; it requires $`\mathrm{Y}`$ to be a $`\mathrm{Z}^0`$ boson and is applied to data taken at all energies. The data used for these searches were recorded by the OPAL detector at centre-of-mass energies ($`E_{\mathrm{cm}}`$) $`88209`$ GeV, the entire energy range achieved at LEP.

These searches are largely motivated by “fermiophobic” scenarios where one of the Higgs bosons decays primarily into a boson pair. In the fermiophobic interpretation, Y would be a $`\mathrm{Z}^0`$ and X a Higgs boson decaying into two photons. Indeed, the Higgs boson predicted in the Standard Model decays into two photons via a quark- or W-boson loop , but with a rate too low for observation of the process at LEP luminosities. Processes $`\mathrm{e}^+\mathrm{e}^{}\mathrm{h}^0\mathrm{Z}^0\gamma \gamma \mathrm{f}\overline{\mathrm{f}}`$ at near-Standard Model production rate and having large di-photon branching ratios have been predicted in a number of alternative theories ; here $`\mathrm{h}^0`$ refers to the lightest neutral boson where extended Higgs sector models are discussed. A particularly natural situation for fermiophobic Higgs bosons occurs in two Higgs doublet models (2HDM)  of “Type-I”, where one Higgs doublet couples only to bosons. Because there are different fermiophobic models, it is not possible to present search results for the entire parameter space of the various theories. In the present paper a benchmark fermiophobic model is defined as having Standard Model production strength and a Higgs boson di-photon branching fraction calculated by turning off the fermion couplings to the Higgs boson in the Standard Model.

The OPAL Collaboration has presented searches similar to those reported here for LEP energies up to $`E_{\mathrm{cm}}`$ = 189 GeV ; this paper extends those searches with the addition of data taken at $`E_{\mathrm{cm}}`$ = $`192209`$ GeV. Fermiophobic Higgs boson searches have also been presented by other LEP collaborations  and by hadron collider experiments . To date, no evidence of a fermiophobic Higgs boson has been seen.

## 2 Data, Simulated Backgrounds and Signals

The data used in this analysis were recorded using the OPAL detector  at LEP. The 1999 data consisted of $`217.0\pm 0.7`$ pb-1 collected at $`E_{\mathrm{cm}}`$ = $`192202`$ GeV. The 2000 LEP data consisted of $`211.1\pm 0.8`$ pb-1 collected at $`E_{\mathrm{cm}}`$ = $`200209`$ GeV, with the majority of the data taken at 205 and 207 GeV. The data sets are summarized in Table 1.

The backgrounds from Standard Model processes were modelled using Monte Carlo simulations at $`\sqrt{s}=192`$, 196, 202, and 206 GeV for the 1999 and 2000 data. Simulated events were processed using the full OPAL detector Monte Carlo  and analysed in the same manner as the data. The full-detector simulations were reweighted for the $`\sqrt{s}`$ distribution of the data using the Monte Carlo generators.

The dominant background to this search arises from the emission of two energetic initial state radiation (ISR) photons. This process was simulated using the KK2f/CEEX  generator with hadronisation and fragmentation by PYTHIA 6.125 . The CEEX modelling of ISR employs full second-order QED corrections to the matrix element, and applies coherent exponentiation of the QED corrections from interference between ISR and final-state radiation. Four-fermion processes were modelled using the grc4f  and KORALW  generators. Two-lepton final states were simulated using BHWIDE , TEEGG  and KORALZ . The NUNUGPV  program was used to generate events of the type $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\gamma \gamma (\gamma )`$. The process $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ was simulated using the RADCOR generator . Tau lepton decays were modelled using Tauola 2.4 .

The process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{h}^0\mathrm{Z}^0`$, $`\mathrm{h}^0\gamma \gamma `$ was simulated for each $`\mathrm{Z}^0`$ decay channel using the HZHA3 generator . For the general search, which is applied to the data taken in 2000 only, the role of the $`\mathrm{Z}^0`$ was replaced by scalar or vector particles having masses from $`10200`$ GeV. Efficiencies for signals were estimated by generating Monte Carlo for both scalar and vector signals in mass steps of 5 GeV; the scalar and vector efficiencies agreed within systematic errors and therefore are not treated as separate cases.

## 3 Event Selection

The analysis described in the following is identical to the one used in the paper for OPAL data taken at 189 GeV . Slightly different analysis cuts were used on the lower energy data sets, as described in the earlier publications .

Events were selected if there were at least two photons recoiling from some other system decaying into one of the following three topologies:

* a $`\mathrm{q}\overline{\mathrm{q}}`$ pair (“Hadronic Channel”), or
* one or two charged leptons (“Leptonic Channel”), or
* a $`\nu \overline{\nu }`$ pair (“Missing Energy Channel”).

Photons were identified as clusters in the electromagnetic calorimeter (EC) which were not associated with tracks if the lateral spread of the clusters satisfied the criteria described in reference . The efficiencies were increased by approximately $`1020`$% by including photon conversions into $`\mathrm{e}^+\mathrm{e}^{}`$ pairs using the methods described in reference .

The dominant background to the searches arises from ISR producing mostly low-energy photons along the beam direction. Therefore, we required the two highest-energy photons in the event to satisfy the following:

* The two photon candidates were required to be in the fiducial region $`|\mathrm{cos}(\theta _\gamma )|<0.875`$, where the polar angle $`\theta _\gamma `$ is the angle of the photon with respect to the $`\mathrm{e}^{}`$ beam direction.
* The highest-energy photon was required to have $`E_{\gamma 1}/E_{\mathrm{beam}}>0.10`$ and the second-highest-energy photon was required to have $`E_{\gamma 2}/E_{\mathrm{beam}}>0.05`$.
* The sum of track momenta and extra electromagnetic cluster energies in a 15 degree cone about the photons had to be less than 2 GeV.

The remaining cuts depend on the particular recoil topologies. In order to assess the background modelling, the photon cuts are not applied until after preselection cuts for the three final state topologies. For all topologies, tracks and EC clusters that are not associated to tracks are required to satisfy the criteria defined in reference . The criteria for the definition of tracks and EC clusters in the OPAL detector are described in reference .

### 3.1 Hadronic Channel

The hadronic channel is characterised by two photons recoiling against a hadronic system. Candidate events were required to satisfy the following criteria:

* The standard OPAL hadronic event preselection in Ref. ; $`R_{\mathrm{vis}}>0.5`$; $`|\mathrm{\Sigma }p_\mathrm{z}^{\mathrm{vis}}|<0.6E_{\mathrm{beam}}`$; and at least two electromagnetic clusters with $`E/E_{\mathrm{beam}}>0.05`$. The quantities $`E_{\mathrm{vis}}`$ and $`\stackrel{}{p}_{\mathrm{vis}}`$ are the scalar and vector sums of track momenta, unassociated EC and unassociated hadron calorimeter cluster energies, and $`R_{\mathrm{vis}}\frac{E_{\mathrm{vis}}}{E_{\mathrm{cm}}}`$. The visible momentum along the beam direction, obtained from the sum of all tracks and unassociated clusters, is denoted by $`|\mathrm{\Sigma }p_\mathrm{z}^{\mathrm{vis}}|`$.
* The photon pair criteria G1$``$G3.
* Photon isolation: both photon candidates were required to satisfy $`p_{\mathrm{T},\mathrm{jet}\gamma }>5`$ GeV, where $`p_{\mathrm{T},\mathrm{jet}\gamma }`$ is the photon momentum transverse to the axis of the closest jet out of two jets formed with the Durham  scheme (excluding the photon pair).
* Photon energy balance: $`(E_{\gamma 1}E_{\gamma 2})/E_\mathrm{o}<0.5`$, where $`E_\mathrm{o}(sM_\mathrm{Z}^2)/(2\sqrt{s})`$ would be the energy of a single photon recoiling from the $`\mathrm{Z}^0`$. This cut discriminates against ISR photon pairs.
* The recoil mass from the di-photon system, $`M_{\mathrm{recoil}}`$, is required to be consistent with the $`\mathrm{Z}^0`$: $`|M_{\mathrm{recoil}}M_\mathrm{Z}|<20`$ GeV (not used in the general search mode).

### 3.2 Charged Lepton Channel

This channel searches for events in the $`\gamma \gamma \mathrm{}^+\mathrm{}^{}`$ final state. Events having only one well-identified lepton are accepted to avoid efficiency loss for lepton tracks at low polar angles. The lepton tracks are treated as jets to include tau lepton final state topologies. Leptonic channel candidates were required to satisfy the following criteria:

* The standard OPAL low multiplicity preselection of Ref. ; $`R_{\mathrm{vis}}>0.2`$; $`|\mathrm{\Sigma }p_\mathrm{z}^{\mathrm{vis}}|<0.8E_{\mathrm{beam}}`$; number of EC clusters not associated with tracks $`10`$; number of tracks $`N_\mathrm{T}`$ satisfies $`1N_\mathrm{T}7`$; at least two electromagnetic clusters with $`E/E_{\mathrm{beam}}>0.05`$.
* The photon pair criteria G1$``$G3.
* For events having only one track and a converted photon, the EC cluster associated with the track must not also be associated with the conversion.
* For events having two or more tracks, the event is forced to have two jets within the Durham scheme, excluding the identified di-photon candidate, and both jets are required to have energies above 3 GeV.
* $`|M_{\mathrm{recoil}}M_\mathrm{Z}|<20`$ GeV (not used in the general search mode).

### 3.3 Missing Energy Channel

The missing energy channel is characterised by two photons and no other significant detector activity. Candidates in the missing energy channel were required to satisfy the following criteria:

* The standard OPAL low multiplicity preselection of Ref. ; the vetoes in Ref.  against cosmic ray and beam-wall/beam-gas backgrounds; number of EC clusters not associated with tracks $`4`$; number of tracks $`3`$; $`|\mathrm{\Sigma }p_\mathrm{z}^{\mathrm{vis}}|<0.8E_{\mathrm{beam}}`$; and at least two electromagnetic clusters with $`E/E_{\mathrm{beam}}>0.05`$.
* The photon pair criteria G1$``$G3.
* $`p_T(\gamma \gamma )>0.05E_{\mathrm{beam}}`$ where $`p_T(\gamma \gamma )`$ is the transverse momentum of the di-photon system; the angle between the two photons in the plane transverse to the beam axis: $`|\varphi _{\gamma \gamma }180^{}|>2.5^{}`$; the polar angle of the momentum of the di-photon system: $`|\mathrm{cos}\theta _{\gamma \gamma }|<0.966`$.
* No track candidates other than those associated with an identified photon conversion.
* Veto on unassociated calorimeter energy: the energy observed in the EC not associated with the two photons is required to be less than 3 GeV.
* $`|M_{\mathrm{recoil}}M_\mathrm{Z}|<20`$ GeV (not used in the general search mode).

## 4 Results

For the 1999 and 2000 data, the numbers of events passing the cuts are listed for the three recoil topologies in Table 2. There are no events in which more than one photon pair satisfying the cuts was found. The numbers of candidates passing cuts are generally in good agreement with the expected numbers of Standard Model backgrounds; this was also the case in earlier OPAL publications for the lower $`E_{\mathrm{cm}}`$ . The one noteworthy discrepancy is for cut C1 in the missing energy channel. In this channel there is a large background from Bhabha electrons lost in the beampipe. The ISR photons for this background have a steeply rising population in the forward direction. Cut C1 is made before the cut on polar angle, and therefore a discrepancy arises because of the steep angular distribution and the inadequate modelling of material in the very low polar angle regions.

Combining both the 1999 and 2000 data in the three topologies, 112 candidates pass the general cuts compared to 118.3$`\pm `$7.9 expected background, and 42 candidates pass the $`\mathrm{h}^0`$$`\mathrm{Z}^0`$cuts compared to 51.9$`\pm `$2.9 expected background.

### 4.1 Systematic Errors

The uncertainty on the modelling of ISR is the most important component of the systematic error because of the irreducible background arising from this process. This uncertainty is estimated from the comparison of data with the Standard Model background simulation for events passing cuts A2, B2, or C2. The shapes of the distributions for $`E_{\gamma 1}`$ and $`E_{\gamma 2}`$ are modelled well by the simulations. The simulations also reproduce well the number of events observed in the three channels combined. For the 1999 and 2000 runs the simulations predict 3.7% and 5.0% fewer events than observed, respectively. The statistical error on the 1999 data is 4%, and similarly for the 2000 data. Modelling of photon conversions has an uncertainty of approximately 1%. Uncertainties on the integrated luminosities of the data sets are negligible compared to the other uncertainties. Combined, these error sources result in a total background uncertainty estimate of 10%. The experimental results which follow are not very sensitive to this number.

The dominant systematic uncertainty for the signal acceptances arises from the photon detection efficiency, primarily due to the simulation of the photon isolation criterion G3 , and is estimated to be 3%. The uncertainty from Monte Carlo statistics is typically better than 4%. A systematic error on the photon energy scale is estimated by comparing the fitted single-photon ISR energy peak to the expected value based on the precisely known beam energy and $`\mathrm{Z}^0`$ mass. For this study a sample of single-photon events was generated and compared to the data for photon energies above 5 GeV and polar angles greater than 25 degrees. This leads to a systematic uncertainty on the di-photon mass of 0.35 GeV at a mass of 100 GeV. The resolution on the di-photon mass ranges from approximately 0.5 GeV at $`m_{\gamma \gamma }`$ = 10 GeV to 2.3 GeV at $`m_{\gamma \gamma }`$ = 100 GeV.

### 4.2 General Search Results

Figure 1 shows the di-photon mass versus the recoil mass for all candidate events passing the general search cuts for the year 2000 data only (where all the data were taken at $`E_{\mathrm{cm}}`$ near 206 GeV). The events at recoil masses near zero are expected from $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ background. This plot also shows no unexpected structure for the lower $`E_{\mathrm{cm}}`$ data. In the absence of an indication for signal, limits are placed on the production at $`E_{\mathrm{cm}}`$$``$206 GeV of a massive state decaying into photon pairs.

For the general search, the system recoiling from the di-photon system is not assumed to be a $`\mathrm{Z}^0`$ and hence the branching fractions $`\mathrm{X}\gamma \gamma `$ and Y into topology A, B or C are not uniquely predicted. Here X is a scalar particle and $`\mathrm{Y}`$ is a scalar or vector particle. Furthermore, X and Y must be a short-lived particles so that they decay near the interaction point. In order to be independent of models we do not combine data from different $`E_{\mathrm{cm}}`$ and therefore we restrict this part of the analysis to the highest energy data, in the $`E_{\mathrm{cm}}`$ range of $`205207`$ GeV; this represents 200.0 pb-1 of the 2000 data. We choose to present upper limits on $`\sigma (\mathrm{e}^+\mathrm{e}^{}\mathrm{XY})\times B(\mathrm{X}\gamma \gamma )\times B(\mathrm{Y}\mathrm{f}\overline{\mathrm{f}})`$ as a function of $`M_\mathrm{X}`$. When presenting production upper limits as functions of $`M_\mathrm{X}`$, we show the limit obtained for the value of $`M_\mathrm{Y}`$ that gives the smallest efficiency in the region $`M_\mathrm{Z}M_\mathrm{X}<M_\mathrm{Y}<\mathrm{E}_{\mathrm{cm}}M_\mathrm{X}`$. The lower bounds on $`M_\mathrm{Y}`$ are used because searches for di-photon resonances at LEP1  have already set good limits on the lower-mass phase space. $`M_\mathrm{X}`$ and $`M_\mathrm{Y}`$ are also required to be above 10 GeV and below 200 GeV in order to allow the decay products to have sufficient energies and momenta to give reasonable search acceptances at $`E_{\mathrm{cm}}`$ = 206 GeV. For a scalar/vector hypothesis for X/Y, the efficiency is found to be the same to within 5% as that for a scalar/scalar hypothesis; the lower of these efficiencies is used in setting the limits. For the lepton search channel, the efficiency for Y $`\tau ^+\tau ^{}`$ is used, as it turns out to have the lowest of the dilepton efficiencies.

The event candidates from $`E_{\mathrm{cm}}`$ = $`205207`$ GeV in the general search are used to calculate 95% CL upper limits on the number of events in 1 GeV \[$`M_\mathrm{X},M_\mathrm{Y}`$\] mass bins. The acceptances used at each 1 GeV mass bin are obtained by interpolation using a 4th-order polynomial fit to the acceptances simulated on a 5 GeV grid. For each \[$`M_\mathrm{X},M_\mathrm{Y}`$\] bin, the 95% CL upper limit on the number of signal events is computed using the frequentist method of reference , which takes into account the predicted Standard Model background. This statistical procedure also incorporates the di-photon mass resolution (typically less than 2 GeV for $`m_{\gamma \gamma }`$$`<`$100 GeV); the limit procedure is valid for resonance states narrower than this resolution. The effect of the 10% systematic error for background modelling is incorporated in the statistical procedure, as is the 4% uncertainty on signal. For these general search limits, an additional systematic uncertainty of 5% is added to the signal uncertainty to account for interpolation error in the efficiency grid (especially near kinematic limits) and for the differences in the acceptance calculations for the scalar versus vector nature of particle Y.

Figure 2 shows the 95% CL upper limits on $`\sigma (\mathrm{e}^+\mathrm{e}^{}\mathrm{XY})\times B(\mathrm{X}\gamma \gamma )\times B(\mathrm{Y}\mathrm{f}\overline{\mathrm{f}})`$. These results are valid independent of the nature of Y, provided it decays to two jets, a lepton pair, or missing energy, and has a width less than or equal to the experimental resolution. Limits of $`2560`$ fb are obtained over $`10<M_\mathrm{X}<180`$ GeV. The limits for the leptonic final state are stronger, except in the case Y couples exclusively to $`\tau ^+\tau ^{}`$ (the final state with lowest acceptance).

### 4.3 Limits on $`\mathrm{h}^0`$$`\mathrm{Z}^0`$ with $`\mathrm{h}^0`$$`\gamma \gamma `$

The distribution of di-photon masses for the $`\mathrm{h}^0`$$`\mathrm{Z}^0`$ search candidates for the 1999 data ($`E_{\mathrm{cm}}`$ = $`192202`$ GeV) and the 2000 data ($`E_{\mathrm{cm}}`$ = $`200209`$ GeV) is shown in Figure 3a together with the simulation of Standard Model backgrounds. The observation of 42 events is in reasonable agreement with the expected background of 51.9$`\pm `$2.9 events. Because the $`\mathrm{h}^0`$$`\mathrm{Z}^0`$ process has a production rate and branching fractions described by theory, the data taken at all LEP energies can be combined in this analysis. Figure 3b shows the distribution of $`m_{\gamma \gamma }`$ for $`E_{\mathrm{cm}}`$ = $`88209`$ GeV. This plot is restricted to $`m_{\gamma \gamma }`$ larger than 20 GeV because there is no background estimate for the low-$`m_{\gamma \gamma }`$ LEP1 data. The figure has no indication of a resonance, and the total of 124 candidates agrees with the predicted background of 135.2$`\pm `$10.8 events.

Also indicated on Figure 3 is the hypothetical signal of a 100 GeV Higgs boson produced at Standard Model strength and decaying into di-photons with a branching fraction of 18% – the fraction predicted in the “benchmark fermiophobic model” calculated by simply turning off the Higgs-fermion coupling. In reality, the fermiophobic Higgs photon branching ratio depends on parameters and details of fermiophobic 2HDM models , so this benchmark is simply a guide to the broad interpretation of the data. Here we use the HDECAY  package to calculate the modified photonic branching fractions.

The events passing all $`\mathrm{h}^0`$$`\mathrm{Z}^0`$ cuts are used to set an upper limit on the di-photon branching ratio for a particle produced in association with a $`\mathrm{Z}^0`$ and having the Standard Model Higgs boson production rate. As described in the previous section, the frequentist method of reference  is used to determine the 95% confidence level upper limit on possible signal events at each di-photon mass. Figure 4 shows the 95% CL upper limit for the di-photon branching ratio obtained by combining the candidate events in 1999 and 2000 data described in this paper with those from OPAL searches at $`\sqrt{s}=88189`$ GeV , where the Standard Model $`\mathrm{h}^0`$$`\mathrm{Z}^0`$ production cross-section is assumed at each centre-of-mass energy. Higgs bosons produced at Standard Model rate and decaying exclusively to di-photons are ruled out at the 95% confidence level over the mass range $`20117`$ GeV. Figure 4 also shows the $`\mathrm{h}^0\gamma \gamma `$ branching ratio computed using HDECAY with the fermionic couplings switched off; the photonic branching fraction falls as the $`\mathrm{W}^+\mathrm{W}^{}`$ and $`\mathrm{Z}^0\mathrm{Z}^0`$ channels become kinematically favourable. A 95% CL lower mass limit for the benchmark fermiophobic Higgs bosons is set at 105.5 GeV, where the predicted branching ratio crosses the upper-limit curve. The median limit one would expect to obtain in an ensemble of experiments in the absence of a signal is 106.4 GeV. The benchmark fermiophobic branching ratios can also be calculated using the HZHA3  generator. HZHA3 produces slightly higher di-photon branching fractions than does HDECAY. The lower mass limit on the benchmark fermiophobic Higgs boson calculated with HZHA3 is 106.3 GeV.

## 5 Conclusions

A search for the production of Higgs bosons and other new particles of width no larger than the experimental resolution and decaying to photon pairs has been performed using $`\mathrm{e}^+\mathrm{e}^{}`$ annihilation data with $`E_{\mathrm{cm}}`$ = $`192209`$ GeV combined with $`88189`$ GeV data from previous OPAL searches. Model independent upper limits are obtained for $`E_{\mathrm{cm}}`$$``$206 GeV on $`\sigma (\mathrm{e}^+\mathrm{e}^{}\mathrm{XY})\times B(\mathrm{X}\gamma \gamma )\times B(\mathrm{Y}\mathrm{f}\overline{\mathrm{f}})`$, where limits of $`2560`$ fb are obtained over $`10<M_\mathrm{X}<180`$ GeV, for $`10<M_\mathrm{Y}<200`$ GeV and $`M_\mathrm{X}+M_\mathrm{Y}>M_\mathrm{Z}`$. The limits are valid for Y either a scalar or vector particle, provided that the Y decays to a fermion pair (interpreted as two jets, a lepton pair, or missing energy).

Using OPAL data from all LEP centre-of-mass energies, model-specific limits are placed on $`B`$($`\mathrm{h}^0\gamma \gamma `$) up to a Higgs boson mass of 117 GeV, provided the Higgs particle is produced via $`\mathrm{e}^+\mathrm{e}^{}\mathrm{h}^0\mathrm{Z}^0`$ at the Standard Model rate. A lower mass bound of 105.5 GeV is set at the 95% confidence level for benchmark fermiophobic Higgs bosons. Similar lower mass limits on benchmark fermiophobic Higgs bosons have been obtained by the other LEP experiments .

Acknowledgements:

We thank A. G. Akeroyd, L. Brücher, and R. Santos for helpful discussions. We particularly wish to thank the SL Division for the efficient operation of the LEP accelerator at all energies and for their close cooperation with our experimental group. In addition to the support staff at our own institutions we are pleased to acknowledge the
Department of Energy, USA,
National Science Foundation, USA,
Particle Physics and Astronomy Research Council, UK,
Natural Sciences and Engineering Research Council, Canada,
Israel Science Foundation, administered by the Israel Academy of Science and Humanities,
Benoziyo Center for High Energy Physics,
Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) and a grant under the MEXT International Science Research Program,
Japanese Society for the Promotion of Science (JSPS),
German Israeli Bi-national Science Foundation (GIF),
Bundesministerium für Bildung und Forschung, Germany,
National Research Council of Canada,
Hungarian Foundation for Scientific Research, OTKA T-029328, and T-038240,
Fund for Scientific Research, Flanders, F.W.O.-Vlaanderen, Belgium.


text:
# The winds of Luminous Blue Variables and the Mass of AG Car

## 1. Introduction

The strong winds of LBVs show a wide variety of mass-loss behaviour. During their S Dor-type variations they expand in radius at approximately constant luminosity. In some cases the mass loss increases while the star expands (e.g. R 71), whereas for others (e.g. R 110) the behaviour is the exact opposite: as the star expands, its mass-loss rate drops. Recent radiation-driven wind models of OBA supergiants show that stars change their wind characteristics at spectral types B1 and A0, where $`\dot{M}`$ jumps upwards by factors of five, due to Fe recombinations. In this poster, we investigate whether these “bi-stability jumps” can also explain $`\dot{M}(T_{\mathrm{eff}})`$ of LBVs.

## 2. Mass Loss Predictions

Our method is outlined in the poster by de Koter & Vink (these proceedings). The full study is presented in Vink & de Koter (2002). Typical LBV results are shown in the first figure, for three values of $`v_{\mathrm{}}`$/$`v_{\mathrm{esc}}`$. It shows that there are ranges in $`T_{\mathrm{eff}}`$ where $`\dot{M}`$ is predicted to increase, and ranges where $`\dot{M}`$ decreases. The decreases are due to a growing mismatch between the positions of the driving lines (mostly in the UV) and the location of the flux maximum which shifts towards the optical for cooler stars. The increases are due to recombinations of Fe iv to Fe iii, and Fe iii to Fe ii. UV observations have shown that these “bi-stability” jumps occur at spectral types B1 and A0 (Lamers et al. 1995).

## 3. Comparison with the observed Mass Loss behaviour of AG Car

The second figure shows a comparison between our models (dotted line) and the H$`\alpha `$ mass loss rates for AG Car (Stahl et al. 2001). It shows that the observed and predicted mass loss agree within 0.1 dex. Note that we have applied a corrective shift of $`\mathrm{\Delta }T_{\mathrm{eff}}=\mathrm{6\hspace{0.17em}000}`$ K to our predictions to account for an inaccurate calculation of the ionisation balance of Fe. The applied shift is consistent with constraints set by observations of supergiants which show that the jumps indeed occur at spectral types B1 and A0. As $`\dot{M}`$($`T_{\mathrm{eff}}`$) displays such a complex behaviour with fluctuations in $`\dot{M}`$ of over more than 0.5 dex, this is a surprisingly good result. It shows that the mass-loss variability of AG Car is due to changes in the ionisation balance of iron.

## 4. The Mass of AG Car

The predicted mass-loss rates are not only a function of $`T_{\mathrm{eff}}`$, but also of Mass. This sensitivity of $`\dot{M}`$ to Mass offers the opportunity to constrain LBV Masses. Comparing the AG Car data with $`\dot{M}`$ predictions for different masses results in a best fit to the present-day Mass of 35 $`M_{}`$ for AG Car.

## 5. Conclusions

* LBV winds are driven by radiation pressure.
* The mass loss behaviour of LBVs (of up to over 0.5 dex) during their S Dor-type variation cycles can be explained by the ionisation and recombination of Fe iv/iii and Fe iii/ii. noindent
* The $`\dot{M}`$($`T_{\mathrm{eff}}`$) behaviour of AG Car can be matched when we adopt a mass of 35 $`M_{}`$.

## References

Lamers H.J.G.L.M.,Snow T.P., & Lindholm D.M., 1995, ApJ 455,269

Stahl O., Jankovics I., Kovacs J., et al., 2001, A&A 375, 54

Vink J.S., & de Koter A., 2002, A&A in press, astro-ph/0207170