text:
# Chandra Observations of the Nucleus of M33

## 1 Introduction

The nearby galaxy M33 ($`d795`$ kpc; van den Bergh 1991) hosts the most luminous steady X-ray source in the Local Group. M33 X-8 dominates the X-ray emission from the galaxy with a 1-10 keV luminosity of about 1.2$`10^{39}`$ erg s-1, ten times more than that of the other X-ray sources in the field. The ultra-luminous X-ray source (ULX) has been detected at this level since its discovery by Einstein. (Long et al., 1981; Markert & Rallis, 1983; Gottwald et al., 1987; Peres et al., 1989; Takano et al., 1994; Long et al., 1996; Parmar et al., 2001). Observations obtained with the ROSAT High Resolution Imager (HRI) showed that X-8 is point-like and coincident with the optical nucleus of M33 at the HRI 5″ resolution (Schulman & Bregman, 1995). These also revealed a 106 day $``$20% modulation of the X-ray flux from M33 X-8 (Dubus et al., 1997).

The nature of the X-ray emitting object is rather puzzling. The measured velocity dispersion of the nucleus places an upper limit of 1500 M on the mass of a putative massive central black hole (Gebhardt et al., 2001). A low luminosity active galactic nucleus (AGN) is therefore unlikely. Indeed, the UV/optical spectrum of the nucleus shows no signs of AGN-type activity and is best modelled by a combination of two starbursts at 40 Myrs and 1 Gyrs (Long, Charles & Dubus 2002 and references therein).

The 106 day modulation suggests a single object is responsible for most of the X-ray emission, perhaps a persistent X-ray binary with a 10 M black hole radiating at the Eddington luminosity. Hubble Space Telescope observations place a limit of 18% on the contribution of a point source to the total far-UV flux from the nucleus (Dubus et al., 1999b). The implied ratio $`L_{\mathrm{opt}}`$/$`L_\mathrm{X}0.05`$ would be in line with expectations from a low-mass X-ray binary. The counterpart would then be mostly lost in the glare of the nucleus. However, there remains a distinct possibility that M33 X-8 is associated with a UV bright star located 1″ to the NNW (Dubus et al., 1999b). The star has a O9III spectrum with no outstanding feature linking it to the powerful X-ray source but cannot be ruled out at the ROSAT HRI resolution (Long, Charles & Dubus, 2002).

Most studies have assumed M33 X-8 is associated with the nucleus and used this to register the X-ray images. Chandra offers the opportunity to locate M33 X-8 to unprecedented accuracy and settle whether this source is within the nucleus or not. Three observations of M33 are available in the Chandra archive (Table 1). We used the first two to accurately position M33 X-8 with respect to HST Planetary Camera (PC) images of the nucleus (§2). The third, non piled-up observation is used to study the X-ray spectrum and timing behaviour of M33 X-8 (§3 and §4). Our results are discussed in the last section.

## 2 Astrometry

The 1-2″ precision of the Chandra pipeline aspect solution111http://cxc.harvard.edu/cal/ASPECT/celmon/index.html is insufficient for our purposes. To improve on this, we identified X-ray sources within the field-of-view and cross-correlated these to optical, IR and radio sources with accurately known positions.

### 2.1 X-ray position of the nucleus

M33 X-8 is severely affected by pile-up in the two 50 ks ACIS observations. The count rate in the ACIS-S observation is 0.5 counts/frame compared to 3.0 counts/frame in the 10 ks non piled-up 10 ks ACIS-S observation (see §3 below). The pile-up fraction is therefore about 84% which makes it impossible to study the radial extension of the source to higher accuracy than the ROSAT HRI. Observations with the Chandra HRC would be able resolve any contribution from extended emission or faint point sources within 5″. The inner pixels of the ACIS-S (Fig. 1) are very piled-up hence have a dearth of counts so that the 0.5-2 keV image is doughnut shaped (standard G02346 events). Going to higher energies reduces the effect and there is a very well defined peak in both ACIS-S and I.

CELLDETECT (CIAO v2.2.0.1222http://cxc.harvard.edu/ciao/) locates M33 X-8 at (4098.1,4057.2)$`\pm `$0.1 pixels in the ACIS-I and at (4063.0,4026.3)$`\pm `$0.5 pixels in the ACIS-S. To check these values, we redefined the center as the location which minimizes the distances of all the events within a 300 pixel wide box roughly centered on M33 X-8 (this assumes circular symetry and avoids the problem posed by a pile-up doughnut). The values obtained by selecting only the events with energies above 1, 2 or 5 keV are all consistent with the CELLDETECT locations and errors which we therefore adopt.

### 2.2 Optical position of the nucleus

The optical position of the nucleus we adopt is the 2MASS position which is tied to the Tycho reference frame. The (compact) nucleus is detected with $`J=12.08\pm 0.04`$ in this survey. At $`\alpha (\mathrm{J2000})1^\mathrm{h}33^\mathrm{m}50\stackrel{s}{.}89`$ and $`\delta 30^{}39\mathrm{}36\stackrel{}{.}75`$ with an uncertainty of 0$`\stackrel{}{.}`$15 this is, as far as we know, the best available absolute position. The position of the close O9III star was then found relative to the nucleus using the UV HST PC image described in Dubus et al. (1999b). We find $`\alpha 1^\mathrm{h}33^\mathrm{m}50\stackrel{s}{.}85`$ and $`\delta 30^{}39\mathrm{}37\stackrel{}{.}65`$ to the 2MASS uncertainty.

As an independent check, we cross-correlated the brightest stars in the various HST PC images of the nucleus (Dubus et al., 1999b) with the M33 star catalogue produced by the DIRECT project (Macri et al., 2001a). The DIRECT positions are tied to the USNO-A2 reference frame with an accuracy of 0$`\stackrel{}{.}`$4-0$`\stackrel{}{.}`$7. We found 4-6 stars (depending upon the sky orientation of the PC chip) that correlated well with DIRECT positions, all lying at the same offset within 0$`\stackrel{}{.}`$1-0$`\stackrel{}{.}`$4 of each other. We found a mean corrected position for the nucleus of $`\alpha 1^\mathrm{h}33^\mathrm{m}50\stackrel{s}{.}90`$ and $`\delta 30^{}39\mathrm{}37\stackrel{}{.}14`$. This is consistent with the 2MASS position within the limits of the DIRECT astrometric solution.

### 2.3 X-ray source list and cross-correlation

We applied CELLDETECT to identify X-ray sources in the two 50 ks observations. The detection treshold was 4$`\sigma `$, resulting in 51 (77) possible sources in the ACIS-I (ACIS-S) dataset. We compared this list with the ROSAT X-ray catalogue of Haberl & Pietsch (2001) to identify previously known X-ray sources. ROSAT and Chandra sources lying within 5″ of each other were considered identical and are thereafter identified by their number in Table 1 of Haberl & Pietsch (2001).

We cross-correlated the X-ray source list with the Tycho-2 (Hog et al., 2000), USNO-A2 (Monet et al., 1998) and 2MASS (Skrutskie et al., 1997) catalogues. We also used a subset (Macri et al., 2001b; Mochejska et al., 2001) of the DIRECT catalogue containing variable optical sources (excluding cepheids) and the radio source catalogue of Gordon et al. (1999). The latter contains optically confirmed supernova remnants (SNR) in M33 and stellar-like background radio sources. The radio positions are accurate with respect to the ICRS frame to 7″/2$`\sigma `$, where 7″ is the beam size and $`\sigma `$ is the signal-to-noise ratio of the source.

We started by searching for all possible counterparts in the catalogues within 10″ of each X-ray source and plotting their pixel offsets relative to the Chandra position. A given X-ray source can be associated with several sources from the same or different catalogues. In the plot, true counterparts will cluster at the offset of the X-ray astrometry (the orientation error and ACIS chip distortion are negligible), while spurious identifications will appear uniformly distributed (the field of view is much greater than the searched correlation radius). The spread of the cluster is a measure of the accuracy of the astrometric solution. There was an obvious concentration of matches around the offset between X-8 and the 2MASS position of the nucleus (Table 2 and Fig. 2). Note that some of the matches may refer to the same object such as the USNO and Tycho match for HP87. None of the DIRECT variable sources matched an X-ray source within 3″ although we note that a star with a 3.4 days optical period (D33J013334.8+303211.4, Mochejska et al. 2001) lies within 10″ of the 3.45 day eclipsing X-ray source X-7 (Dubus et al., 1999a).

### 2.4 Position of the nuclear X-ray source

The clustering of matches around the offset between the 2MASS position of the nucleus and the Chandra position of X-8 is particularly clear for the radio correlations (Fig. 2). For uniformly distributed sources, the probability $`p`$ of having an X-ray source within $`d`$ arcsec of a possible counterpart is $`p(d/D)^2`$ if $`dD`$ and where $`D`$ is the radius of the field of regard. Correlating $`N`$ X-ray sources with $`M`$ possible counterparts in the field of view, the probability $`P`$ of observing $`C`$ or more sources within $`d`$ arcsec of a counterpart is

$$P=1\underset{i=0}{\overset{C1}{}}\left(\genfrac{}{}{0pt}{}{MN}{i}\right)p^i(1p)^{MNi}$$

For $`C=8`$ matches with the radio catalogue, $`d=3`$″, $`N=32`$, $`M=126`$ and $`D8`$′, the probability of seeing a cluster of 8 associations within 3″ is $`<\mathrm{\hspace{0.33em}10}^{11}`$. The same calculation for the ACIS-S observation ($`C=7`$, $`N=57`$, $`M=117`$) yields a probability $`<\mathrm{\hspace{0.33em}10}^8`$. We conclude the radio associations are very likely to be true counterparts (in contrast with the USNO, 2MASS and Tycho associations for which the same calculation does not rule out chance superpositions).

We now estimate the probability that the offsets given by the radio counterparts are consistent with the offset between the 2MASS position of the nucleus and X-8 ($``$ in Fig. 2) or the nearby star position and X-8, plotted as a $``$ in Fig. 2). A $`\chi ^2`$ test on the measured distances between the radio and nucleus offsets yields a probability 0.28 (0.53) that these are consistent for the ACIS-I (ACIS-S) data. The probability for the radio offsets to be consistent with the nearby star are $`210^6`$ (ACIS-I) and $`410^4`$ (ACIS-S). The ACIS-S result is less conclusive than the ACIS-I because of the higher uncertainty in the position of X-8 in the ACIS-S ($`\pm `$0.5 pixel compared to $`\pm `$0.1 pixel, see §2.1). We conclude the nearby star is ruled out as a possible counterpart at the 4.6$`\sigma `$ level and that M33 X-8 is at the 2MASS position of the nucleus with an uncertainty of 0$`\stackrel{}{.}`$6 (the mean weighted distance of the radio offsets for both ACIS-I and S).

There is a radio source associated with the nucleus. This source (GDK102) has $`F_{20\mathrm{c}\mathrm{m}}=0.6\pm 0.1`$ mJy, $`F_{6\mathrm{c}\mathrm{m}}=0.2\pm 0.1`$ mJy is located 0$`\stackrel{}{.}`$6 away from the 2MASS position (and 0$`\stackrel{}{.}`$4 away from the mean radio offset) and is a possible counterpart to M33 X-8. Using the results of the above analysis, we note the association of X-ray source HP84 with 2MASS 0133418+303849 ($`J=12.62\pm 0.02`$) is likely to be real. HP84 is classified as a supersoft X-ray source by Haberl & Pietsch (2001).

## 3 X-ray Spectrum

We performed spectral analysis using data from the ACIS-S observation ObsId 787. We extracted 27259 counts from a circle within 20″ of the source position for a total of 3.06 counts/frame. The source was 7.7′ off-axis, at which the the 50% enclosed energy function (EEF) at 1.5 keV is 4″ vs. 0.3″on-axis, resulting in the source covering an area larger by a factor $``$170, to 0.02 counts/frame/PSF, at which the estimated pileup fraction should be negligible ($`<`$1%). We confirmed this by comparing the distribution of photon event grades (for grades 0, 2, 3, 4 and 6) below 2 keV and in the 2-8 keV range with those of a non-piled up X-ray observation (of ZW3146, ObsId 1651, Seq. Num 800119) and found that they were roughly consistent, though not identical; however, both were similar to the branching ratio observed from a non piled-up observation of PG1634-706 333http://asc.harvard.edu/cal/Hrma/hrma/misc/oac/psf2/index.html, and so pileup does not appear to be significant during this observation. Background counts were taken from an annulus of inner- and outer- radii of 25″ and 100″. The expected number of background counts in the source region is 167$`\pm `$3. The background subtracted source countrate is then 2.92+/-0.02 c/s.

Using XSPEC v11.0.1 (Arnaud, 1996), we produced a pulse-height-analyser (PHA) spectrum of the extracted counts data between 0.6-9 keV, binning at 130 eV between 0.6 and 2.0 keV, at 260 eV between 2.0 and 5.4 keV, and with two energy bins between 5.4 and 9.0 keV. We included a 5% systematic uncertainty in the fits. Best-fit parameters for acceptable models are in Table3.

An absorbed powerlaw spectrum (wabs\*powerlaw) is rejected, with a value of reduced chi-squared statistic $`\chi _\nu ^2`$=4.12 for 23 degrees of freedom (dof), with a null-hypothesis probability of $`10^{10}`$. The power-law spectrum over-predicts counts in the 2-4 keV range, and under-predicts them in the 4-9 keV range. An absorbed blackbody spectrum is also statistically rejected ($`\chi _\nu ^2=21.4`$ for 23 dof).

We find a moderately acceptable fit ($`\chi _\nu ^2`$=1.35) for an absorbed thermal bremsstrahlung spectrum. An absorbed disc blackbody spectrum produces a better fit to the data ($`\chi _\nu ^2`$=0.95), with an X-ray column density $`N_H=(6\pm 1)\times 10^{20}\mathrm{cm}^2`$, consistent with galactic absorption only 444W3NH at http://heasarc.gsfc.nasa.gov/ (Dickey & Lockman, 1990). The disc blackbody inner temperature $`kT_{\mathrm{in}}`$ and inner-radius $`R_{\mathrm{in}}\sqrt{\mathrm{cos}\theta }`$ are identical to the values found in earlier observations with ASCA (Takano et al., 1994) and BeppoSAX (Parmar et al., 2001). The flux (corrected for absorption) of the best-fit spectrum is 2.0$`\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$, for a luminosity 1.5$`\times 10^{39}\mathrm{erg}\mathrm{s}^1`$ (0.5-10 keV), consistent (including spectral uncertainty) with previously measured values. Unlike previous work, a power-law component is not required by the data. When we include a power-law component in addition to the absorbed disc-blackbody model, holding $`\alpha =2`$ fixed, the 90% confidence upper-limit for the flux of this component is 1.1$`\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$.

## 4 X-ray Variability

We examined the data for long-term noise-like time variability. We binned the 27259 counts detected in the full Chandra pass band, with 1.04104 sec time resolution, and took a Fourier transform of the data. With the resulting Fourier components, we obtained a power-density spectrum (PDS; see Lewin et al. 1988 for the use of PDS for broad-band timing analysis) which we rebinned logarithmically. Examination of the resulting PDS shows no evidence of power above that expected from Poisson noise. Modelling the power as a power-law of slope=1 plus a constant (held fixed at the Poisson level) in frequency ($`P(\nu )=A\nu ^1+C`$), we find an upper-limit to the root-mean-square variability of $`<`$3% (90% confidence; 0.0001-0.5 Hz).

When we binned the data by a factor of 1, 10 and 100, (in which the average counts per bin are 2.938, 29.4, and 294 counts), we find that the maximum counts per bin was 11, 48 and 330 counts, which is consistent with the maximum number of counts expected from Poisson fluctuations. Thus, we find no evidence for intensity variability on timescales between 1 and 10,000 sec.

## 5 Discussion

The high luminosity of M33 X-8 is not sufficient to rule out a neutron star X-ray binary since magnetospheric accretion can result in super-Eddington luminosities. This was observed e.g. in the high mass X-ray binaries (HMXBs) SMC X-1 or A0535-66. A HMXB would have been an attractive possibility should M33 X-8 have been associated with the nearby 09III star. The Chandra localisation of X-8 rules this out (§2). X-8 could still be a HMXB hiding in the nucleus but this would be surprising since the compactness and spectra of the nucleus make it unlikely that it hosts $`>10M_{}`$ stars (Long, Charles & Dubus, 2002); a HMXB evolves on the uncomfortably short thermal timescale of the donor star ($`<10^5`$ yrs for $`M_2>10M_{}`$); the X-ray spectrum of X-8 is unlike that of typical HMXBs which have flat power laws in the 0.5-10 keV range; and the power spectrum shows no evidence for pulsations.

On the other hand, the X-ray spectrum and absence of short-term X-ray variability below 1 Hz are highly reminiscent of the high or very high states of galactic black hole candidates (Tanaka & Lewin, 1995). In both states the spectrum is dominated by a disc blackbody with a weak contribution from a steep power law. The very high state shows 3-10 Hz QPOs with a few % rms superposed on 1-10% rms flat top noise in the 1-10 Hz range. There is little variability in the high state. GRS 1915+105, a microquasar in our Galaxy with a $`14`$ M black hole accreting from a $`1.2`$ M K giant star (Greiner, Cuby & McCaughrean, 2001), also has such soft, low variability states (Belloni et al., 2000). The optical emission from such a system would be lost in the nucleus of M33. Since its discovery in 1992, GRS 1915+105 has been steadily detected with X-ray luminosities around or in excess of $`10^{39}`$ erg s-1.

The analogy with GRS 1915+105 is even more attractive considering the 0.6 mJy (respectively 0.2 mJy) radio source detected by the VLA at 20 cm (6 cm; Gordon et al., 1999) which we show to be consistent with the ultra-luminous X-ray source (§2.4). GRS 1915+105 has been known to show steady radio emission around 30-50 mJy at 2 cm, framed by powerful ejection episodes where the emission from the relativistic jets becomes $``$ 500 mJy (Fender et al., 1999). The radio source in M33 is more powerful than that; but a small change in the beaming angle of the jet ($`70^{}`$ from the line-of-sight for GRS 1915+105) could easily increase the observed radio flux. Follow-up observations of the radio source are needed to test this hypothesis.

## ACKNOWLEDGMENTS

We gratefully acknowledge use of the Vizier catalogue (CDS, France), of data products from the 2MASS survey, and of archival HST and Chandra observations obtained from the STScI and SAO websites.


format: text

text:
# Paraconsistency of Interactive Computation Originally published in proc. PCL 2002, a FLoC workshop; eds. Hendrik Decker, Dina Goldin, Jørgen Villadsen, Toshiharu Waragai (http://floc02.diku.dk/PCL/).

## 1 Introduction

The paradigm shift from algorithms to interactive computation captures the technology shift from mainframes to workstations and now to wireless devices and intelligent appliances, from number-crunching to embedded systems and graphical user interfaces, and from procedure-oriented to object-based and distributed computation. Interactive computation, modeled by interactive systems, is a more powerful paradigm of problem solving than algorithmic computation, modeled by Turing Machines.

The view that TMs completely express the intuitive notion of computing is a common misinterpretation of Church-Turing Thesis. This view has led computer science theorists to assume that questions of expressiveness of finite computing agents had been settled once and for all, so that exploration of alternative models of computability has seemed unnecessary. The claim that interactive computing agents are more expressive than Turing Machines, meaning able to solve a larger set of computational problems, requires a reexamination of fundamental assumptions about models of computation, about computation itself, and about computational logics.

In section 2, we introduce interactive computation and demonstrate that it is more powerful than algorithmic computation of Turing Machines, in that it can solve a wider class of computational problems. In section 3, we show that classical first-order logic cannot model interactive computation. In section 4, we examine areas of computing where the principle of consistency plays a major role; in each case, these areas are fundamentally not interactive. Finally, in section 5, we consider embedding logic models in the broader context of interactive models, and conclude that paraconsistency is necessary to model interaction.

## 2 The power of interactive computation

In this section, we introduce interactive computation and demonstrate that it is more powerful than the algorithmic computation of the Turing Machines, in that it can solve a wider class of computational problems.

### 2.1 Interactive computation

The notion of algorithms well precedes computer science, having been a concern of mathematicians for centuries. It can be dated back to a 9th century treatise by a Muslim mathematician Al-Koarizmi, after whom algorithms were named.

> Algorithm: systematic procedure that produces – in a finite number of steps – the answer to a question or the solution of a problem .

Algorithmic computation, modeled by Turing Machines, was adopted as the traditional model of computation when computer science emerged as a new academic discipline in the 1960s.

> Algorithmic computation: the computation is performed in a closed-box fashion, transforming a finite input, determined by the start of the computation, to a finite output, available at the end of the computation, in a finite amount of time.

Whereas algorithmic computation is like the data processing mode associated with punch cards of the 1950s, today’s computational systems are interactive. In these systems, which are expected to remain up and running continuously, input is consumed and output is generated throughout the computation. More importantly, there is an interdependence of input and output: output can influence later input, and vice versa.

### 2.2 The *WH* problem

Let us consider an automatic car whose task is to drive us across town from point *W* (work) to point *H* (home); we shall refer to it as the WH problem. The output for this problem should be a time-series plot of signals to the car’s controls that enable it to perform this task autonomously. At issue is what form the inputs should take.

The car can be equipped with a map of the city. In the algorithmic scenario, where all inputs are provided a priori of computation, this map needs to include every pothole and even grain of sand along the road. By the principles of chaotic behavior, such elements can greatly affect the car’s eventual course – like the Japanese butterfly that causes a tsunami on the other end of the world.

In a static world, such a map is in principle obtainable; but ours is a dynamic environment. The presence of mutable physical elements such as the weather affect the car’s course both directly (e.g., the wind blowing at the car) and indirectly (e.g. sand shifting location in the path of the car as a result of blowing wind). It is arguable whether the elements of the weather can be precomputed to an accuracy required to produce the proper answer to the *WH* problem.

We can remain optimistic until we remember that the world also includes humans, as pedestrians or drivers. To avoid collisions, we must precompute the motion of everyone who might be near the car as it drives. An assumption that human actions can be computed ahead of time is tantamount to an assertion of fatalism – a doctrine that events are fixed in advance so that human beings are powerless to change them – clearly beyond the purview of any computer scientist. Therefore, we must assume it cannot be done.

> The claim that computational tasks situated in a dynamic world that includes human agents are solvable algorithmically is tantamount to a assertion of fatalism.

### 2.3 The power of interaction

Nevertheless, the *WH* problem is solvable – interactively. In this scenario, the inputs, or percepts , consist of a stream of images produced by a video camera mounted on the car, as it is driving from *W* to *H*. The signals to the car’s controls are generated on-line in response to these images, to avoid steering off the road or running into obstacles. This change in the scenario is akin to taking the blindfolds off the car’s electronic driver, who was driving from memory and bound to a precomputed sequence of actions. Now, he is aware of his environment and uses it to guide his steering.

Note that a change of approach from algorithmic to interactive such as above involves much more than a restructuring of inputs. The nature of computation itself is at stake. Algorithmic computation is off-line; it takes place before the driving begins. Interactive computation is on-line; it takes place as the car drives. The values of inputs and outputs for interactive computation are interdependent; decoupling them, such as by replacing the videocamera with a prerecorded videotape of the road, will usually lead to an incorrect solution.

This interpretation of the concept of computation is different from the one that Church and Turing had in mind for their famous thesis, where they specifically referred to computation of functions over finite strings that encode natural numbers (i.e., algorithmic computation). Interactive computation falls outside the bounds of the Church Turing thesis. The *WH* example shows that interaction pulls us out of the Turing tarpit, by providing solutions to problems that are not solvable in the algorithmic setting.

> By allowing us to solve computational tasks that cannot be solved algorithmically, interaction is shown to be a more powerful computational paradigm.

## 3 Logic and interactive computation

In this section, we discuss why classical first-order logic cannot model interactive computation.

### 3.1 Tradeoffs between reasoning and modeling

Logic comprises both syntactic rules for noninteractively (in a closed-box fashion) proving theorems from axioms by rules of inference (proof theory) and semantic notions of soundness and completeness that relate syntactic processes of proof to semantic properties of a modeled domain in which formulae are interpreted (model theory). A logic is sound if all its theorems are true, and complete if all true assertions of the modeled domain are theorems.

In sound and complete formal systems such as first-order logic, every true assertion about the domain has a finite proof that is constructed by applying the rules to the axioms. With the number of rules and axioms also finite, the resulting set of assertions must be recursively enumerable (RE).

> Necessary condition for completeness: A sound and complete first-order logic (SCL) can model only domains with a countable set of properties.

Completeness is a desirable mathematical property that certifies formalizability of the modeled domain but by the same token restricts expressiveness. We may wish to model classes of systems (mathematical, computational, or physical) which have a non-RE set of properties. Modeling is observational rather than constructive, based on greatest fixpoints rather than least fixpoints . The gap between least- and greatest-fixpoint semantics is also the gap between operational (algorithm) and denotational (observation) semantics. It is also the same as the gap between deduction and abduction.

Greatest fixpoints allow us to define larger domains. Lacking a constructive foundation, we cannot systematically enumerate the members of such domains. However, it can be observed whenever two members are distinct; distinguishability certificates are always finite. The bisimulation relation for labeled transition systems is an example of a relation with greatest fixpoint semantics and finite distinguishability certificates.

> Sufficient condition for incompleteness: Any system with an uncountable set of properties cannot be expressed by a sound and complete logic.

As a result, there are tradeoffs between reasoning and modeling: purely rule-based syntactic reasoning is too weak to capture all semantic properties of mathematical, computational, or physical systems. Model builders must choose between complete formalizability and modeling power. Turing Machines focus on formalizability at the expense of modeling power.

(As an aside, this focus on formalizability is intentional; TMs were invented as part of Turing’s proof that Hilbert’s challenge to formalize mathematics cannot be carried out. The halting problem presented a constructive counterexample to this challenge . Turing himself never claimed that TMs formalize all of computation, and in fact made claims that suggested the opposite .)

### 3.2 Incompleteness of interactive systems

Godel showed that arithmetic over the integers cannot be completely formalized by logic, because not all its properties are syntactically expressible as theorems. We extend Godel’s reasoning to show that interactive and empirical systems are likewise incomplete because they have too many properties to be expressible as theorems of a sound and complete logic.

As discussed above, sound and complete first-order logics (SCLs) have an RE set of theorems and can formalize only semantic domains with a countable number of distinct properties. Incompleteness occurs when the number of true facts or observable properties of a system cannot be recursively enumerated by theorems. For example, arithmetic over the integers is incomplete because its true properties cannot be recursively enumerated; that is, for any sound formalization we can find a true property of the integers that cannot be proved.

For interactive systems or interaction machines (IMs), properties correspond to observations, and the number of properties corresponds to the number of distinct observations. Algorithms and TMs have only a countable set of properties, while IMs have an open-ended set of properties that is generally not countable and therefore cannot be expressed as theorems of a complete first-order logic. Godel’s incompleteness result holds for any domain whose true properties cannot be recursively enumerated, including empirical domains modeled by IMs.

Expressiveness for IMs is defined in terms of the ability of observers to make observational distinctions . Observational expressiveness allows behavior of algorithms and IMs to be measured by the same metric. The class of systems or programs which admits finer observational distinctions is the one capable of a greater range of behaviors, and hence able to express solutions to a larger set of problems. From this viewpoint, TM behavior is completely characterized by a single observation, while IM behavior may depend on unbounded sequences or patterns of observations. This gives another proof of greater expressive power of interaction.

Completeness restricts the expressiveness of domains modeled by an SCL: domains whose properties are not countable are necessarily incomplete. IMs, based on infinite streams of inputs and outputs created in partnership with an (uncomputable) environment, have an uncountable set of computations: the distinct streams of an IM are not enumerable.

By viewing Godel’s incompleteness result as a corollary of the more general result that SCLs cannot model systems with an uncountable set of properties we gain insights into the nature of incompleteness. Incompleteness is seen to be a ubiquitous phenomenon that applies not only to mathematical systems like the integers but also to IMs.

> The set of behaviors of an interactive system cannot in general be formalized by any sound and complete logic.

Incompleteness shows the limitations of syntactic formalisms in expressing semantic behavior. The completeness/incompleteness dichotomy distinguishes algorithms from IMs, closed from open systems, and rationalist from empiricist models . Formal incompleteness of computing systems is related to descriptive incompleteness of physical models of the real world.

Interactive problem solving (section 2.2) is a first-class form of computation that should be included in any intuitive notion of computability. Hence, logics need to go beyond soundness and completeness to capture computation.

### 3.3 Logic for physical theories

The relatively recent development of paraconsistent logic challenges the logical principle that contradictory premises lead to an explosion of meaningless conclusions.

> Paraconsistency: Let $``$ be a relation of logical consequence. Let us say that $``$ is explosive iff for every formula $`A`$ and $`B`$, $`\{A,\neg A\}B`$. Classical logic and most other standard logics are explosive. A logic is said to be paraconsistent iff its relation of logical consequence is not explosive .

The most telling reason for paraconsistent logic is the fact that there are non-trivial theories which are inconsistent. Some of the best examples of such theories come from the history of physics.

The interactive view of computing allows us to bridge the gap between closed-box algorithmic computing and real physical systems. Interactive computation spans the range from closed-box algorithmic computation at one end to situated environment-aware computing agents (e.g. the car in the *WH* problem) in the middle, to the behavior of real world objects at the other extremum.

The notion of observation as a metric of expressiveness and equivalence for IMs is a natural extension of the role of observation in the formulation of physical theories . The idea that physical objects are not completely describable or knowable but that they may have describable parts or views is a basic tenet of the scientific method. An analogous idea for IMs is that there is no silver bullet – when testing a software system, one can never be sure of having found all the bugs.

> Just as physical theories provide an important motivation for paraconsistent logics, so does the theory of interactive computation.

## 4 Consistency in computation

In this section, we examine some of the areas of computing which have been founded on the principle of consistency, i.e., that it is impossible for $`A`$ and $`\neg A`$ to coexist. We show that in each case, these areas are fundamentally not interactive.

### 4.1 Logic programming

Logic programming corresponds to closed-system computing, whereas interactive systems are inherently open. The Japanese 5th generation computing project that aimed to reduce computation to logic programming proved a failure. It was shown in  that this was due not to lack of cleverness on the part of logic programming researchers, but to the theoretical impossibility of such a reduction.

The key argument is the inherent incompatibility between reactiveness of interactive computation, realized by committed choice at each interaction step, and logical completeness, realized by backtracking. Committed choice is inherently incomplete because commitment cuts off branches of the proof tree that might contain the solution. However, commitment is essential to interactive computation, since execution of an observable action by a system cannot be reversed. The power of interaction can be achieved only by sacrificing logical completeness. Pure logic programming is inherently too weak to serve as a model for interactive computation.

### 4.2 Data Management

Despite significant overlaps, the research areas of databases (DB) and information systems (IS) are based on fundamentally different assumptions, due to the different jobs they perform. The “job of a database” is to store data and answer queries. This entails addressing issues like data models, schema design, handling distributed data, maintaining data consistency, query evaluation, etc.

Given a query and a database instance, the behavior of a database system consists of answering queries, which is algorithmic in nature. In particular, a database management system implements an algorithmic mapping of the form:

> (DB contents, user query) $``$ (DB contents’, query answer).

On the other hand, the job of an IS is to provide information-based services, which entails considerations that span the life cycle of the larger system. Services over time are interactive in nature, involving user sessions with one or more users.

It is not surprising then that consistency is taken for granted in database research, given its algorithmic nature. However, most large contemporary database systems are to some degree IS’s, providing additional services beyond simply storing data and running queries and updates . When lines between the two are blurred in this fashion, the tension between consistency and interaction can result.

### 4.3 Artificial Intelligence

Classical logic has always played a prominent role in Artificial Intelligence (AI), where a typical problem would be to form a plan or to solve a problem, given a complete description of the “world” in question. Closed-box computation was the modus operandi, even in problems like rearranging blocks in Shakey’s world  that were meant to offer a simplified version of situated real-world behavior.

Many researchers came to realize that algorithmic approaches cannot scale up to solutions for tasks involving the real world , like the *WH* problem; a paradigm shift was needed. The new agent-based approach to AI  is interactive. The change in approach has reinvigorated this field, breaking the inertia that gripped it in the ’80s, and restoring our optimism in the promises that AI was making when originally founded.

Not surprisingly, classical logic does not play the same prominent role in the new AI. The emphasis has shifted to probablistic approaches, with such models as neural networks, hidden Markov models and belief networks . We expect that paraconsistent logics can also play a significant role in the new AI.

## 5 Extending logic to interaction

By embedding logic models in the broader context of interactive models, we gain insight into the limits of logic in modeling autonomous external worlds. We also understand why paraconsistency is necessary to model interaction.

The number of interpretations of function and predicate symbols of first-order logic (fol) is nonenumerable. However, the set of theorems for a fixed interpretation is recursively enumerable (section 3.1). By assuming a static interpretation of the world that is fixed prior to the start of the computation, algorithmic models have a recursively enumerable set of behaviors. By contrast, the nonenumerability of interpretations of fols reflects the nonenumerable external worlds (environments) of interactive computation, which are neither static nor known a priori. Interactive models replace complete semantics of fols by weaker abstractions that sacrifice completeness for expressiveness.

> Conjecture: Late (dynamic) binding of interpretations in logic is a form of interactive modeling.

Interaction extends models so that the interpretation of nonlogical symbols evolves during inference to take account of new data. However, interactive discovery of facts negates the fundamental property that true facts always remain true. Nonmonotonic logics permit reasoning about incompletely described modeled worlds. As the system gains information that contradicts earlier inferences (beliefs), it may retract these inferences.

The contrast between the constructive paradigm of algorithmic models and the observational paradigm of interaction is the contrast between deduction and abduction (section 3.1). Abduction is open-ended and subject to revision; abductive reasoning is paraconsistent in that it is possible for incompatible inferences to be drawn from the same premises.

Abduction and nonmonotonicity are not only desirable, but required for interactive computation, where input, inference, and action are interleaved. Paraconsistency, the embodying of contradictory information, is necessary for abductive and nonmonotonic reasoning, and hence for interaction.

## 6 Conclusion

The goal of computational logic is to allow us to model computation as well as to reason about it. A useful computational logic must be able to model interactive computation. We showed that first-order logic cannot model interactive computation due to its incompleteness. Interactive computation is necessarily paraconsistent, able to model both a fact and its negation, due to the role of the world (environment) in determining the course of the computation. We noted that interactive computation shares much in common with scientific theories, providing a common motivation for the development of paraconsistent logic. We concluded that paraconsistency is a necessary property for a logic that can model interactive computation.

## Acknowledgements

We thank David Keil for his careful reading of earlier versions of this work.


text:
# Using Balmer line profiles to investigate convection in A and F stars

## 1. Introduction

Balmer lines are an important diagnostic of stellar atmospheric structure since they are formed at a wide range of depths within the atmosphere. The depth of formation of H$`\alpha `$ is higher than that of H$`\beta `$, thus observations of these profiles provide useful diagnostics (e.g. Gardiner 2000). Balmer profiles are relatively insensitive to surface gravity for stars cooler than $``$8000 K (Gray 1992, see also Heiter et al. 2002), whilst sensitive to the treatment of atmospheric convection (e.g. van’t Veer & Mégessier 1996, Castelli et al. 1997, Gardiner 2000, Heiter et al. 2002). For stars hotter than $``$8000 K the profiles are sensitive to both effective temperature and surface gravity. However, provided we know surface gravity from some other means (e.g. from eclipsing binary systems), we can use them to determine effective temperature.

In previous work, Smalley & Kupka (1997) found no significant systematic problems with $`uvby`$ and fundamental (and standard) stars. In fact, $`uvby`$ was found to be very good for obtaining $`T_{\mathrm{e}ff}`$ and $`\mathrm{log}g`$. Using H$`\alpha `$ and H$`\beta `$ profiles, Gardiner et al. (1999) found that both the Canuto & Mazzitteli (1991, 1992) and standard Kurucz (1993) mixing-length theory without overshooting (see Castelli et al. 1997) are both in agreement to within the uncertainties of the fundamental stars. Overshooting models were always clearly discrepant. However, Gardiner et al. (1999) found some evidence for significant disagreement between all treatments of convection and fundamental values around 8000$``$9000 K. In this region the effects of $`\mathrm{log}g`$ cannot be ignored, and the majority of the $`T_{\mathrm{e}ff}`$ stars did not have fundamental values of $`\mathrm{log}g`$. We have used binary systems with fundamental values of $`\mathrm{log}g`$, determined fundamental values of $`T_{\mathrm{e}ff}`$ and compared the results with those from fitting models to Balmer-line profiles.

## 2. Effective temperatures of binary systems

Eclipsing binary systems provide ideal test stars for comparing to models, since they enable us to obtain fundamental values of $`T_{\mathrm{e}ff}`$ and $`\mathrm{log}g`$. We can obtain fundamental values of $`T_{\mathrm{e}ff}`$ provided we know the apparent angular diameter and total integrated (bolometric) flux at the Earth. In the case of binary systems, where there are no direct measurements of angular diameters, we can obtain them from the stellar radius and the parallax of the system (Smalley & Dworetsky 1995).

Available spectrophotometry was taken from various sources. Unfortunately, not all the systems have enough high-resolution spectrophotometry. In these cases, however, we have at least $`UBV`$ or $`UBVRI`$ colours, which were used to estimate the optical fluxes. Near-infrared fluxes were taken from the Gezari et al. (1987), 2MASS and DENIS catalogues. The lack of available fluxes for binary systems is something that the Citadel ASTRA Spectrophotometer (Adelman et al. 2002) could address.

In most cases the $`T_{\mathrm{e}ff}`$ values obtained using the Hipparcos parallaxes are in agreement with that obtained from other methods (e.g. Infrared Flux Method, $`uvby`$ photometry, etc.) and other determinations (e.g. Jordi et al. 1997, Ribas et al. 1998). However, three systems were anomalous and these are discussed in detail in Smalley et al. (2002).

## 3. Observations

The H$`\alpha `$ and H$`\beta `$ observations were made at the Observatorio del Roque de los Muchachos, La Palma using the Richardson-Brealey Spectrograph on the 1.0m Jacobus Kapteyn Telescope (JKT) in October/November 1997. A 2400 l $`mm^1`$ holographic grating was used and a 1124 $`\times `$ 1124 pixel Tek CCD, giving a resolution of 0.4Å fwhm. Further H$`\alpha `$ and H$`\beta `$ observations were made at the Mount Stromlo Observatory, Australia in February 2000 using the Cassegrain Spectrograph on the ANU 74 inch telescope. A 1200 l $`mm^1`$ blazed grating was used, giving a resolution of 0.35Å fwhm. Additional H$`\beta `$ profiles were taken from the work of Smalley & Dworetsky (1995).

The data reduction of the profiles taken in 1997 and 2000 was performed using the Starlink echomop software package. In most cases the final spectra had a signal-to-noise ratio in excess of 100:1. Instrumental sensitivity variations were removed from the H$`\alpha `$ profiles by comparing to observations of stars with intrinsically narrow Balmer profiles, for example early-B or O type stars and G type stars, and the H$`\beta `$ profiles were normalized such that the observed profile of Vega agreed to a model with $`T_{\mathrm{e}ff}`$=9550 K, $`\mathrm{log}g`$=3.95, \[M/H\]=$``$0.5 (Castelli & Kurucz 1994) and the standard profiles of Peterson (1969).

## 4. Effective temperatures from Balmer line profiles

The observed Balmer line profiles are fitted here to model spectra to compare the derived $`T_{\mathrm{e}ff}`$ with that from fundamental methods. The following convection models were used, using solar-metallicity Kurucz atlas models:

Standard atlas9 (Kurucz 1993) models using mixing length theory (MLT) without convective overshooting. The value of the MLT parameter $`\alpha `$ is the standard value of 1.25.

Standard atlas9 models using MLT without convective overshooting. The value of the MLT parameter $`\alpha `$ is 0.5.

Standard atlas9 models using MLT with approximate overshooting. The value of the MLT parameter $`\alpha `$ used is 1.25.

Standard atlas9 models using MLT with approximate convective overshooting. The value of the MLT parameter $`\alpha `$ used is 0.5.

Modified atlas9 models using the Canuto & Mazzitelli (1991,1992) model of turbulent convection.

The synthetic spectra were calculated using uclsyn (Smalley et al. 2001) which includes Balmer line profiles calculated using VCS Stark broadening and metal absorption lines from the Kurucz & Bell (1995) linelist. This routine is based on the balmer routine (Peterson 1969). The synthetic spectra were normalized $`\pm 100`$ Å to match the observations. The values of $`T_{\mathrm{e}ff}`$ were obtained by fitting model profiles to the observations using the least-square differences.

Figures 1 & 2 show the variation of $`\mathrm{\Delta }`$$`T_{\mathrm{e}ff}`$ = $`T_{\mathrm{e}ff}`$(Balmer) $``$ $`T_{\mathrm{e}ff}`$(fund) against $`T_{\mathrm{e}ff}`$(fund) for H$`\alpha `$ and H$`\beta `$, respectively, for the 5 convection models listed above. To within the uncertainties, the CM results show no significant variation with $`T_{\mathrm{e}ff}`$(fund) for either H$`\alpha `$ or H$`\beta `$. The discrepancy around 8000 K noted by Gardiner et al. (1999) is not evident. Even the two anomalous H$`\alpha `$ points just hotter than 8000 K, for V624 Her, can be brought into agreement if the IRFM $`T_{\mathrm{e}ff}`$ is used (Smalley et al. 2002). The MLT\_noOV results are in broad agreement with those for CM, but with the $`\alpha `$=0.5 models giving better agreement around 8000 K relative to $`\alpha `$=1.25 and CM models. Contrary to Gardiner et al. (1999), who reported that F-type stars might require models with $`\alpha `$1.25 (see their Fig. 9), we find that the binary systems do not support this. Overall, $`\alpha `$=0.5 models are preferred to those with higher values. The MLT\_OV models are generally more discrepant, yielding too high values of $`T_{\mathrm{e}ff}`$ (and even larger ones for H$`\beta `$, if $`\alpha `$=1.25 is used rather than 0.5), as found previously by Gardiner et al. (1999). Note also the systematic difference between H$`\alpha `$ and H$`\beta `$ for $`\alpha `$=1.25 MLT\_noOV models, which is even more pronounced for the MLT\_OV models.

## 5. The apparent A-star anomaly

The use of stars with fundamental values of both $`T_{\mathrm{e}ff}`$ and $`\mathrm{log}g`$ has failed to support the apparent anomaly around 8000 K found by Gardiner et al. (1999). However, there were too few stars within the $`T_{\mathrm{e}ff}`$ range 8000–9000 K to fully explore this region. Gardiner et al. (1999) found that four fundamental $`T_{\mathrm{e}ff}`$ stars also showed the anomaly: $`\beta `$ Leo, $`\alpha `$ Oph, $`\alpha `$ Aql, $`\alpha `$ PsA. In order to be sure that there is no anomaly in the Balmer line profiles, we need to explain why these stars might appear anomalous.

Table 2 summarizes the values of $`T_{\mathrm{e}ff}`$ obtained from CM $`uvby`$ photometry, the IRFM and by fitting to H$`\alpha `$ and H$`\beta `$ profiles. We have allowed both $`T_{\mathrm{e}ff}`$ and $`\mathrm{log}g`$ to vary in order to obtained the best least-squares fit (see Figures 3 & 4). Values of $`\mathrm{log}g`$ are also given as obtained from $`uvby`$ photometry. We have also included $`\gamma `$ Gem which is just hotter than 9000 K, but the results are in agreement with its fundamental and IRFM $`T_{\mathrm{e}ff}`$ values.

The rapidly rotating star $`\alpha `$ Aql has recently been studied by van Belle et al. (2001) using interferometry. Their analysis revealed the oblateness of the star and a new determination of fundamental $`T_{\mathrm{e}ff}`$ = 7680$`\pm `$90 K. This is significantly cooler than the previous determination, but in accord with that inferred from the IRFM. As such, the $`T_{\mathrm{e}ff}`$ from H$`\alpha `$ and H$`\beta `$ are no longer significantly discrepant. It is certainly possible that revision to the other fundamental stars could occur once new interferometric measurements are obtained, especially $`\alpha `$ Oph which has a similar $`v\mathrm{sin}i`$ and might be expected to exhibit significant oblateness. Thus, the anomalies for these two stars can be explained in terms of their rapid rotation.

The two other stars, $`\beta `$ Leo and $`\alpha `$ PsA, have lower $`v\mathrm{sin}i`$ values, but are two most discrepant stars in the Gardiner et al. (1999) sample. Unless the fundamental values are truly wrong there must be some other reason for the discrepancy. The IRFM values both point to a slightly cooler $`T_{\mathrm{e}ff}`$, but even then the discrepancy is $``$300 K. However, in this temperature region the Balmer lines are near their maximum strength and sensitive to $`\mathrm{log}g`$. It is certainly possible that a small error in adopted $`\mathrm{log}g`$ could lead to a large error in $`T_{\mathrm{e}ff}`$ obtained from Balmer profiles. In addition, the Balmer profiles change little with relatively large changes in $`T_{\mathrm{e}ff}`$. Thus, we conclude that the two stars are not discrepant, due to the low sensitivity of Balmer lines with respect to changes in $`T_{\mathrm{e}ff}`$ and both sensitivity to, and the uncertainty in, the surface gravity for these stars. However, it must be noted that the $`\mathrm{log}g`$ obtained from the Balmer lines for these stars is systematically lower than that obtained from $`uvby`$ photometry.

In general, for stars hotter than 8000 K the sensitivity to $`\mathrm{log}g`$ prevents us from using them to obtain values of $`T_{\mathrm{e}ff}`$ to the accuracy required for the present task, unless we have accurate fundamental values of $`\mathrm{log}g`$. However, until we do have stars with accurate fundamental $`\mathrm{log}g`$ values, we cannot be totally sure that there is not a problem with the model predictions in this $`T_{\mathrm{e}ff}`$ region.

## 6. Conclusion

Balmer line profiles have been fitted to the fundamental binary systems. To within the errors of the fundamental $`T_{\mathrm{e}ff}`$ values, neither the H$`\alpha `$ or H$`\beta `$ profiles exhibit any significant discrepancies for the CM and MLT without approximate overshooting models. As in previous work, the MLT with overshooting models are found to be discrepant. Moreover, there are no systematic trends, such as offsets, between results from H$`\alpha `$ and H$`\beta `$ as long as $`\alpha `$ in MLT models is chosen small enough (e.g. 0.5). The discrepancies exhibited by the fundamental $`T_{\mathrm{e}ff}`$ stars in Gardiner et al. (1999) can be explained by rapid rotation in two cases and by the fact that the Balmer profiles become sensitive to $`\mathrm{log}g`$ and less sensitive to $`T_{\mathrm{e}ff}`$ in the other two cases. However, for the time being the lack of any stars with fundamental values of both $`T_{\mathrm{e}ff}`$ and $`\mathrm{log}g`$ in this region precludes the conclusion that there is not a problem with the models in the $`T_{\mathrm{e}ff}`$ range 8000 $``$ 9000 K.

Full details of this work are given in Smalley et al. (2002).

### Acknowledgments.

This work has made use of the hardware and software provided at Keele by the PPARC Starlink Project. Friedrich Kupka acknowledges support by the project Turbulent convection models for stars, grant P13936-TEC of the Austrian Fonds zur Förderung der wissenschaftlichen Forschung.

## References

Adelman S.J., Gulliver A.F., Smalley B., Pazder J.S., Younger P.F., Boyd L., Epand D., 2002, Poster presented at this Symposium

Canuto V., Mazzitelli I., 1991, ApJ, 370, 295

Canuto V., Mazzitelli I., 1992, ApJ, 389, 724

Castelli F., Kurucz R.L., 1994, A&A, 281, 817

Castelli F., Gratton R.G., Kurucz R.L., 1997, A&A, 318, 841

Gardiner R.B., Kupka F., Smalley B., 1999, A&A, 347, 876

Gardiner R.B., 2000, Ph.D. Thesis, University of Keele

Gezari D.Y., Schmitz M., Mead J.M., 1987, Catalog of Infrared Observations, NASA Ref. Publ. 1196

Gray D.F., 1992, Obs. and Analysis of Stellar Photospheres, CUP

Heiter U., Kupka F., van’t Veer-Menneret C., Barban C., Weiss W.W., Goupil M.-J., Schmidt W., Katz D., Garrido R., 2002, A&A, (in press)

Jordi C., Ribas I., Torra J., Giménez A., 1997, A&A, 326, 1044

Kurucz R.L., 1993, Kurucz CD-ROM 13: ATLAS9, SAO, Cambridge, USA.

Kurucz R.L., Bell B., 1995, Kurucz CD-ROM 23: Atomic Line List, SAO, Cambridge, USA.

Peterson D.M., 1969, SAO Spec. Rept., No. 293

Ribas I., Giménez A., Torra J., Jordi C., Oblak E., 1998, A&A, 330, 600

Smalley B., Dworetsky M.M., 1995, A&A, 293, 446

Smalley B., Gardiner R.B., Kupka F., Bessell M.S., 2002, A&A, (submitted)

Smalley B., Kupka F., 1997, A&A, 328, 349

Smalley B., Smith K.C., Dworetsky M.M., 2001, UCLSYN Userguide

van Belle G.T., Ciardi D.R., Thompson R.R., Akeson R.L., Lada E.A., 2001, ApJ, 559, 1155

van’t Veer-Menneret C., Megessier C., 1996, A&A, 309, 879

Vidal C.R., Cooper J., Smith E.W., 1973, ApJS, 25, 37


text:
# Distribution of the second virial coefficients of globular proteins

## Abstract

George and Wilson \[Acta. Cryst. D 50, 361 (1994)\] looked at the distribution of values of the second virial coefficient of globular proteins, under the conditions at which they crystallise. They found the values to lie within a fairly narrow range. We have defined a simple model of a generic globular protein. We then generate a set of proteins by picking values for the parameters of the model from a probability distribution. At fixed solubility, this set of proteins is found to have values of the second virial coefficient that fall within a fairly narrow range. The shape of the probability distribution of the second virial coefficient is Gaussian because the second virial coefficient is a sum of contributions from different patches on the protein surface.

PACS: 87.14.Ee, 87.15Nn.

Protein crystallisation is an important problem yet our grasp of the details of how it occurs is very poor. Proteins need to be crystallised from solution in order to determine their structure via X-ray crystallography . The crystallisation presumably starts with heterogeneous nucleation of the crystalline phase in the protein solution, but there has been no systematic experimental study of this, as far as the author is aware. Without an understanding of how proteins crystallise, protein crystallisation is almost totally ad hoc: essentially the only way to know if a protein will crystallise under a certain set of conditions is to try it. It would be enormously useful if we could predict the conditions under which a protein was most likely to crystallise. Here by protein we mean globular protein, which are proteins that are soluble in solution, as opposed to membrane proteins which exist embedded in a surfactant bilayer. The hope that it is possible to predict the conditions that promote crystallisation motivated George and Wilson  to look at the values of the (osmotic) second virial coefficient of a number of proteins under the conditions where they were crystallised. They found that the second virial coefficient was always negative and lay within what they called ‘a fairly narrow range’. If we ignore outliers then second virial coefficients gathered together by Haas and Drenth , and converted to reduced units by Vliegenthart and Lekkerkerker , lie in the range $`8`$ to $`40`$, in units of the volume of the protein; see Table V of Ref. .

The simplest explanation of this range is that the upper limit is set by the requirement that the attractive interactions be strong enough to pull the molecules into a crystal from a dilute solution. The lower limit is set by the dynamics of the solution, if the attractive interactions are too strong the protein molecules tend to aggregate irreversibly and this aggregation preempts and prevents crystallisation. Testing these explanations is all but impossible due to our poor understanding of crystallisation so we turn to a well-defined, and easily calculated, property of a protein solution: its solubility. We consider the solubility of the protein, i.e., the concentration of protein in the solution which coexists with the crystal, in preference to the process of crystallisation. We ask the question: For a given solubility, say 5% by volume, what is the distribution of values that we expect for $`B_2`$? If we have 1000 proteins, say, all with the same solubility, then is their distribution of values of $`B_2`$ very broad, or is it narrow? What is the shape of the distribution, i.e., what is its functional form?

The distribution of values of $`B_2`$ of a large number of proteins defines a probability distribution function $`P(B_2)`$. We will consider a constraint, that of fixed solubility, and so will obtain a probability distribution function that also depends on this constraint. We are inspired to study this function by a range of work on protein solutions and crystals  that has shown that protein-protein interactions are well described by a sum over contacts between the proteins. Where by a contact between two proteins, we mean that a specific patch on the surface of one protein approaches closely, a couple of Å, to a specific patch on the surface of the other protein. Now, if these contacts are more-or-less independent then we expect them to contribute essentially independently to $`B_2`$. But if $`B_2`$ is the sum of many independent contributions for each protein, then we know the form of the distribution function $`P(B_2)`$: it is a Gaussian. The central limit theorem states that the probability distribution function of some property $`Y`$, which is a sum over a large number of independent random variables, is a Gaussian . Also of course the more independent patches there are on the surface the narrower will be the distribution of values of $`B_2`$.

The physical picture is that the surface of a protein has a number of patches on its surface. Under the conditions where the protein’s solubility is low, these patches attract each other. The strength of each patch attraction is then a random variable selected from some distribution. It is a random variable if the strength of the attraction of one patch on the surface is independent of the attraction of any of the other patches. This physical picture is very simple and is of course approximate but if the correlations between the various patches are weak then $`P(B_2)`$ should be approximately Gaussian. Of course this will also apply to the probability distribution function of any other variable which is a sum of more-or-less independent contributions from the surface patches. We do not know whether the rate of crystallisation is such a variable.

Before we consider our model it is worthwhile noting that within biology there is a move away from studying proteins one at a time to studying them en masse, e.g., studying the complete proteome of an organism. Where the proteome is defined as being the complete set of proteins possessed by an organism. This follows on from work on establishing the complete genome of a number of organisms . Although here we consider sets of proteins which are just those we want to crystallise, and so can come from a number of different organisms, in the future the solubilities and virial coefficients of complete proteomes could be considered.

The model is chosen to be as simple as possible, while incorporating the patchy nature of the surface of proteins together with the variability in the interactions from protein to protein. Thus for simplicity we chose a lattice model. The lattice is cubic and each protein occupies 8 lattice sites arranged 2 by 2 by 2, see Fig. 1. We make the model 2 sites across to reduce the range of the attraction, which is 1 site, to half the diameter of the hard core. The model ‘proteins’ can rotate, and so have 24 distinct orientations. Each of the 6 faces of the cube has a patch, labeled $`i=1`$ to 6, with patches 1 to 4 clockwise around a loop of 4 of the faces, and patches 5 and 6 on the remaining 2 faces. The interactions between model proteins are pairwise additive and consist of 2 parts. The first is simply an excluded volume interaction: 2 proteins cannot overlap. The second is that if the faces of 2 proteins are in contact there is an energy of interaction between the 2 touching patches of the 2 proteins. By in contact we mean that all of the 4 lattice sites on the face of one protein are in contact with one of the lattice sites of the face of the other protein, in other words the faces must overlap completely otherwise the energy of interaction is taken to be zero. The touching patches are those on the faces of the 2 proteins that face each other. If the $`i`$th patch of one protein is adjacent to the $`j`$th patch of another protein then there is an interaction energy of $`ϵ_{ij}`$. Different proteins will have the same excluded volume interactions but the set of interactions $`ϵ_{ij}`$ will be different to represent the different surfaces of different proteins. The interactions form a symmetric square matrix, $`ϵ_{ij}=ϵ_{ji}`$.

Thus, a protein is specified by giving values to the 21 distinct $`ϵ_{ij}`$; these are composed of 36 interactions, consisting of 6 like interactions, $`i=j`$, and 30 unlike interactions, $`ij`$, with 15 of these related to the other 15 by symmetry. The energies $`ϵ_{ij}`$ determine the phase behaviour via Boltzmann weights $`b_{ij}=\mathrm{exp}(ϵ_{ij}/kT)`$, where $`k`$ is Boltzmann’s constant and $`T`$ is the temperature. We will in fact deal mainly with these weights not with the energies themselves.

For the values of the $`b_{ij}`$ for the patch-patch interactions we again choose the simplest possible model and neglect any correlations Each of the 21 distinct $`b_{ij}`$ is almost a stochastic variable described by a probability function $`p(b)`$, i.e., the probability that $`b_{ij}`$ lies between $`b`$ and $`b+\mathrm{d}b`$ equals $`p(b)\mathrm{d}b`$. They are almost but not quite uncorrelated because we want the crystalline phase of all our model proteins to be a simple cubic lattice with 1 protein molecule per unit cell. The crystal is close-packed; there is a protein on every lattice site. In such a crystal, 3 different bonds are formed: the $`13`$, $`24`$ and $`56`$ bonds and the energy $`e_x`$ is

$`e_x`$ $`=`$ $`ϵ_{13}+ϵ_{24}+ϵ_{56}`$
$`{\displaystyle \frac{e_x}{kT}}`$ $`=`$ $`\mathrm{ln}\left(b_{13}b_{24}b_{56}\right)`$ (1)
We want this to be the ground state and so having generated the 21 distinct Boltzmann weights, we find the 3 largest of these, call them $`b_{\alpha \beta }`$, $`b_{\gamma \delta }`$ and $`b_{\delta \zeta }`$ and then we interchange $`b_{13}`$ and $`b_{\alpha \beta }`$, $`b_{24}`$ and $`b_{\gamma \delta }`$ and $`b_{56}`$ and $`b_{\delta \zeta }`$. In this way we ensure the simple cubic lattice is the ground state of the model protein. Performing this swapping procedure introduces correlations of course, but they are minor as we have ordered only 3 of the 21 weights.

The entropy of the crystalline phase is zero, as the cubes can neither translate nor rotate, so its free energy is equal to its energy, $`e_x`$. At low temperature, the pressure will be low and so the chemical potential will be closely equal to the free energy, which we have already said is equal to the energy. Thus the chemical potential in the crystalline phase $`\mu _xe_x`$. The chemical potential in a dilute solution, which is well described by an ideal gas, is

$$\mu _{ig}=kT\mathrm{ln}(\rho /24),$$
(2)
where the 24 comes from the rotational entropy $`k\mathrm{ln}24`$ of a freely rotating cube with all 6 faces distinct. $`\rho `$ is the number density, the number of proteins per lattice site. When a dilute solution coexists with the crystal then we can treat the solution as an ideal gas and find the density of this solution by equating $`\mu _{ig}`$ of Eq. (2) with $`\mu _x`$, which is closely equal to $`e_x`$ of Eq. (1). So, the density of the solution that coexists with the crystal, i.e., the solubility of the protein is

$$\rho _s=24\mathrm{exp}(e_x/kT).$$
(3)
A more useful measure of density is the volume fraction $`\varphi =8\rho `$ which the is the fraction of the sites in the fluid occupied by the proteins. Thus we will work in terms of the solubility volume fraction $`\varphi _s=8\rho _s`$.

The second virial coefficient of a continuum model is defined by

$$B_2=\frac{1}{2}d𝐫_{12}d\omega _1d\omega _2\left[\mathrm{exp}\left(u(𝐫_{12},\omega _1,\omega _2)\right)1\right],$$
(4)
where $`𝐫_{12}`$ is the distance between the centres of mass of two molecules, and $`\omega _1`$ and $`\omega _2`$ are the orientations of molecules 1 and 2, respectively. The integrations over the angles are normalised. $`u(𝐫_{12},\omega _1,\omega _2)`$ is the energy of interaction of pair of molecules, as a function of their separation and orientation. For a lattice model the integrals are replaced by sums and we have

$$B_2=\frac{1}{1152}\underset{\alpha =1}{\overset{24}{}}\underset{\beta =1}{\overset{24}{}}\underset{\left\{𝐫_{12}\right\}}{}\left[\mathrm{exp}\left(u(𝐫_{12},\alpha ,\beta )\right)1\right],$$
(5)
where the three sums are, in order, over the 24 orientations of molecule 1, over the 24 orientations of molecule 2 and over the separation of the two molecules. The factor of 1152 comes from the factor of half and two factors of 24 from the normalisation of the sums over orientation.

The second virial coefficient $`B_2`$ of the protein-protein interaction can be calculated once the set of weights $`b_{ij}`$ are specified. The exponential factor in Eq. 5 equals 0 if the two molecules overlap and equals one of the $`b_{ij}`$ when the molecules occupy adjacent sites. Thus we have

$$B_2=\frac{1}{2}\left[27\frac{1}{6}\underset{i=1}{\overset{6}{}}\underset{j=1}{\overset{6}{}}\left(b_{ij}1\right)\right],$$
(6)
where the first term inside the brackets comes from excluded volume interactions and the second from the interactions between touching patches. The number 27 comes from the fact that each model protein excludes other proteins from a cube of 3 by 3 by 3 lattice sites. Thus, in the high temperature limit where $`b_{ij}=1i,j`$, $`B_2=B_{2hc}=27/2`$. The sums over 24 orientations reduce to sums over 6 orientations as rotating either of the 2 molecules around the axis joining their centres does not change the energy. The factor in front of the double sum is a normalisation factor of $`1/36`$ times the 6 possible lattice sites that one molecule can occupy and be adjacent to the other molecule.

Finally, for the purposes of performing example calculations we will take the distribution function $`p(b)`$ to be a top-hat function of mean $`\overline{b}`$ and width $`b_w`$.

$$p(b)=\{\begin{array}{cc}0& b<\overline{b}b_w/2\\ b_w^1& \overline{b}b_w/2<b<\overline{b}+b_w/2\\ 0& b>\overline{b}+b_w/2\end{array}.$$
(7)
With our model defined we can generate a protein by generating values for the 21 distinct Boltzmann weights for the patch-patch interactions, $`b_{ij}`$, according to the probability distribution function $`p(b)`$. Repeating this procedure many times will generate a set of proteins, each protein having a distinct array of interactions, $`b_{ij}`$. This set can then be analysed to look, for example, for correlations between different properties.

Before we look at a set of proteins with solubilities in a fixed, small, range let us look at just a set of proteins, with a range of solubilities. Fixing the mean patch-patch interaction weight $`\overline{b}=9`$ and the width $`b_w=16`$ so that the Boltzmann weights lie in the range $`1`$ to 17, we have generated a set of 10,000 proteins. Humans have about 10 times as many different proteins as this, bacteria typically less than half this number . The mean weight was chosen to give a distribution centred around $`\varphi _s0.05`$. Of course increasing $`\overline{b}`$ will shift the distribution to lower values of $`\varphi _s`$ and decreasing $`\overline{b}`$ will shift the distribution to higher values of $`\varphi _s`$. Their solubilities and second virial coefficients are plotted as a scatter plot in Fig. 2. Clearly, there is a correlation between the second virial coefficient of a protein and its solubility: the more negative is $`B_2`$ the lower is the solubility on average, but there is also considerable scatter.

Now, restricting the solubility of a protein to lie in the narrow range $`\varphi _s=0.05\pm 0.01`$, we have calculated the probability distribution of second virial coefficients and plotted the result in Fig. 3. The curve is almost Gaussian; a Gaussian is a good fit to the data. The standard deviation is close to $`0.2`$. We have used a top-hat probability distribution function for the Boltzmann weights of patch-patch interactions, Eq. 7, but of course our finding of a Gaussian distribution is insensitive to the exact form of this distribution function.

There is relatively little data for proteins where both the second virial coefficient and the solubility have been measured. Extensive results are available for lysozyme , but there are results for both the virial coefficient and the solubility for only a couple of other proteins. Rosenbaum et al.  plot the reduced second virial coefficient against solubility of lysozyme under a number of different conditions, $`\gamma `$-crystallin and BPTI (bovine pancreatic trypsin inhibitor), their Fig. 4. With only 3 proteins no attempt can be made at assessing a probability distribution. However, in the wake of the sequencing of whole genomes there is a drive towards so-called ‘high throughput’ methods which can rapidly assess the properties of large numbers of proteins. If such a method could be developed for second virial coefficients the results could be compared to our predictions. For a given protein there is clear experimental evidence that as conditions are varied so as the make the second virial coefficient more negative the solubility decreases . For our model, this corresponds to increasing the mean Boltzmann weight of the attractions, $`\overline{b}`$, which will decrease both $`B_2`$ and the solubility, whatever the distribution of weights for the patch-patch attractions.

We have swapped 3 of the patch-patch interaction energies to force the ground state to be a simple cubic lattice. Because of this, when the attractions are made strong enough, $`\overline{b}`$ sufficiently positive, our model proteins all have low solubilities. Thus they are presumably representative of proteins that are easily crystallisable. If we had not swapped the 3 interaction energies, and we assume that the only possible crystalline phase is the simple cubic, then some of our model proteins will be highly soluble: they will not crystallise from dilute solutions.

In conclusion, the surfaces of proteins are patchy and mediate short-range interactions, short with respect to the size of the protein, which is a few nms. A quantity such as the second virial coefficient is an integral (sum for lattice models) over the contribution of the core of the protein, over the contributions of any attractions between patches of the surface, and over the contribution of any longer-ranged interactions such as an overall electrostatic repulsion. If the patches are independent or almost independent, then the distribution function for their total contribution to the second virial coefficient, of a large number of proteins of a similar size will tend to be Gaussian. A straightforward consequence of the central limit theorem. The central limit theorem applies to a variable, here the second virial coefficient, that is the sum of a large number of random variables with the same mean and variance. The suggestion is that if experiments are done on a large set of proteins all at the same solubility and all with similar sizes, then this fixed solubility will set a rough scale for the mean strength of the interactions. Then if proteins are highly modular, the predominant variation from protein to protein will come from variations in the sum of the patch-patch Boltzmann weights, and so will have a roughly Gaussian distribution. With a width roughly equal to the the square root of the number of independent patches times the width in the distribution of the Boltzmann weight of a single patch. Thus, a second virial coefficient which is a sum over contributions from a number of independent patches will tend to have a narrow distribution of values. This may be partly responsible for the fairly narrow range of values of the second virial coefficient over which most proteins crystallise  but crystallisation is both complex and poorly understand so there are almost certainly other factors. By highly modular we mean that the patches on a protein’s surface are close to being independent. If they were truly independent we would have $`<b_{ij}b_{i^{}j^{}}>=\overline{b}^2`$, $`iji^{}j^{}`$, where $`<>`$ indicates an average over all the proteins, and the interactions $`ij`$ and $`i^{}j^{}`$ are any two patch-patch interactions of the same protein. We are assuming that violations of this equality are weak.

It is a pleasure to acknowledge discussions with J. Cuesta, D. Frenkel and P. Warren.


text:
# Untitled Document

Quest For Unification

Edward Witten1 Supported in part by NSF Grant PHY-0070928.

Institute For Advanced Study, Princeton NJ 08540 USA

Let us begin by recalling how the standard model of particle physics looked before grand unification was proposed.2 This article is based on my Heinrich Hertz lecture at SUSY 2002 at DESY, June, 2002. Focus, in particular, on the fermions. The left-handed fermions of one generation can be displayed as follows:

$$\left(\begin{array}{c}\nu \\ e^{}\end{array}\right)_1,\left(\begin{array}{c}u\\ d\end{array}\right)_{1/3},\overline{u}_{4/3},\overline{d}_{2/3},e_2^+.$$

Here I have explicitly displayed the pairs $`(\nu ,e^{})`$ and $`(u,d)`$ as weak doublets, but I have not explicitly shown the color quantum numbers of $`u,d`$ and $`\overline{u},\overline{d}`$. The subscript is the hypercharge quantum number $`Y`$.

This structure looks strange at first sight. Part of the strangeness, relative to the way physics was formerly understood, is that the left- and right-handed fermions transform differently under the standard model gauge symmetries. With time, physicists came to see this as a virtue: it means that the quarks and leptons cannot have bare masses and can only gain mass from the Higgs mechanism. So it explains why they cannot be much heavier than the $`W`$ and $`Z`$ bosons (which gain mass from the same Higgs mechanism) and in particular why they are light compared to hypothetical fundamental scales in physics such as the Planck scale of gravity.

But there is another side to the strangeness in the structure of the standard model. A single generation of quarks and leptons is made by putting together quite a few odd bits and pieces, with strange fractions, in particular, for the hypercharge quantum numbers. Anomaly cancellation depends on obscure calculations like

$$\mathrm{Tr}Y^3=2\left(1\right)^3+6\left(1/3\right)^3+3\left(4/3\right)^3+3\left(2/3\right)^3+2^3=0.$$

There is no virtue in this kind of strangeness.

So it was an illumination in 1973 when Georgi and Glashow (building on prior work of Pati and Salam on quark-lepton unification) unveiled the $`SU(5)`$ grand unified theory, or GUT. In this theory, all known elementary particle gauge forces were interpreted as part of a single underlying $`SU(5)`$ gauge force. Indeed, $`SU(5)`$ is the smallest and most obvious simple or unified gauge group in which one can embed the standard model gauge interactions. The embedding is made via a block diagonal ansatz that is based on the fact that $`3+2=5`$:

$$\left(\begin{array}{cc}SU\left(3\right)& \\ & SU\left(2\right)\end{array}\right).$$

(Hypercharge is generated in $`SU(5)`$ by a traceless diagonal matrix that commutes with $`SU(3)\times SU(2)`$.) In this framework, a standard model generation is reduced to two pieces, the $`\mathrm{𝟏𝟎}`$ and the $`\overline{\mathrm{𝟓}}`$, schematically

$$\left(\begin{array}{ccccc}0& \overline{u}& \overline{u}& u& d\\ & 0& \overline{u}& u& d\\ & & 0& u& d\\ & & & 0& e^+\\ & & & & 0\end{array}\right)\left(\begin{array}{c}\overline{d}\\ \overline{d}\\ \overline{d}\\ \nu \\ e^{}\end{array}\right).$$

The fractions are explained as consequences of $`SU(5)`$ group theory, and the verification of anomaly cancellation is greatly abridged. The standard model has never looked the same since the unified $`SU(5)`$ model was proposed.

To achieve these attractive results, it was necessary, as we see in eqn. (1), to place quarks and leptons (and their antiparticles) in the same representations. As in the model of Pati and Salam, this led to a prediction of proton decay. In fact, the model predicts new gauge bosons $`X`$ and $`Y`$, corresponding to the off-diagonal blocks in eqn. (1); they mediate processes such as $`pe^+\pi ^0`$.

There are a few immediate problems. If the strong interactions are unified with the weak and electromagnetic interactions, why are they so much stronger? And will the proton be sufficiently long-lived? These issues were addressed in 1974 in a celebrated paper by Georgi, Quinn, and Weinberg, who calculated the “renormalization group running” of the strong, weak, and electromagnetic couplings. They found that unification was possible, provided that the weak mixing angle had the right value $`\mathrm{sin}^2\theta _W.20`$ and the unification scale was about $`10^{15}`$ GeV. These were spectacular results, since the value of the weak angle was about right, and the value of the unification scale was fortuitous. In fact, $`10^{15}`$ GeV was close enough to the Planck scale to suggest a unification with gravity, and big enough to make the proton long-lived. A proton lifetime of about $`10^{30}`$ years was predicted, large enough to be compatible with experiment but small enough to be observable. The model would have failed if the computed unification scale were significantly less than about $`10^{15}`$ GeV (because the inferred proton lifetime would have been too short) and would have been implausible if the inferred unification scale were much greater than $`10^{19}`$ GeV.

A GUT scale of $`10^{15}`$ GeV may not seem like it is really very close to a Planck scale of $`10^{19}`$ GeV. But this calculation should really be viewed on a log scale, since renormalization group running is logarithmic and what is computed is really the logarithm of the GUT scale. On a log scale, this early computation did give something pretty close to the logarithm of the Planck mass, and it was the first time that any sort of particle physics computation gave a result for any characteristic scale of particle physics phenomena that was at all close to the Planck scale.

It was nice in $`SU(5)`$ to reduce the mess of a standard model generation to just two pieces. But can one do better? It was soon seen (by Georgi and by Fritzsch and Minkowski) that in the larger group $`SO(10)`$, all quarks and leptons of one generation fit neatly in a single irreducible representation. There is a price, though: one has to add a left-handed anti-neutrino. Since it is a standard model singlet, it is natural for it to get a GUT scale mass, and this led to the idea of the “see saw” mechanism for neutrinos (due to Gell-Mann, Ramond, and Slansky and to Yanagida). A model like this will, in the basis $`\left(\begin{array}{c}\nu _L\\ \overline{\nu }_L\end{array}\right)`$, give a mass matrix for the left-handed neutrino and anti-neutrino of the form

$$\left(\begin{array}{cc}0& m\\ m& M\end{array}\right).$$

Here $`m`$ comes from the electroweak Higgs effect, and $`M`$ can be of order the unification scale $`M_{GUT}`$. If we take $`m`$ to be of order $`M_Z`$ and $`M=10^{15}`$ GeV, we get a light neutrino mass in the range $`m_\nu M_Z^2/M(10^2\mathrm{GeV})/10^{15}\mathrm{GeV}=.01`$ eV. This estimate, which was made in the late 1970’s, gave a big impetus to the search for neutrino masses and oscillations. Of course, in making this estimate, it is not clear just what we should put in either the numerator or the denominator, and since experiments of twenty years ago were not sensitive to neutrino masses of $`.01`$ eV, there was considerable interest in variants in which the neutrino masses would be somewhat bigger, for example because $`M`$ arises from loop effects and is smaller than $`M_{GUT}`$.

Having come this far, can we go farther and (i) unify three generations of quarks and leptons in one irreducible representation of a larger gauge group? or (ii) unify Higgs bosons with gauge fields, or with quarks and leptons? The answer to these questions was “no” in the framework of four dimensional grand unification.

For example, there is no candidate grand unified or GUT gauge group that puts several chiral families in one irreducible representation – without “antifamilies” of the opposite chirality. Of course, experimental bounds on such antifamilies have become progressively tighter, and from a theoretical point of view, the possibility that they would combine with the ordinary families and get a large bare mass makes them seem unattractive.

Likewise, in four-dimensional GUT’s, one cannot really unify Higgs particles with quarks and leptons or gauge bosons (or for that matter, unify quarks and leptons with gauge bosons). The closest try uses supersymmetry plus the $`E_6`$ model that I will get to shortly.

The other recognized problem of GUT’s in this period was the “gauge hierarchy problem,” which is the question of why $`SU(2)\times U(1)`$ breaking is so weak compared to GUT symmetry breaking. In the context of the $`SU(5)`$ model, why are $`M_W`$ and $`M_Z`$ so much less than $`M_X`$ and $`M_Y`$?

In this brief review of the GUT theories of the 1970’s, I have stressed the $`SU(5)`$ and $`SO(10)`$ models, which made sense of the fermion quantum numbers and led to predictions of proton decay and neutrino masses. Are there any bigger groups that can teach us more?

Personally, I would say that there is no four-dimensional GUT model that does better – but there is one more model worthy of note. This is the $`E_6`$ model, introduced by Gursey, Ramond, and Sikivie. What is $`E_6`$? In gauge theory, we need to pick a gauge group, which, if we wish to achieve unification, should be a simple gauge group, and moreover should be compact (so that all gauge bosons have positive kinetic energy). There is a nice classification of these groups. First, there are three infinite families, $`SU(N)`$, $`SO(N)`$, and $`Sp(N)`$. Two groups from these infinite families – $`SU(5)`$ and $`SO(10)`$ – are used in the models that I have mentioned so far. Describing nature by a group taken from an infinite family does raise an obvious question – why this group and not another? In addition to the three infinite families, there are five exceptional Lie groups, namely $`G_2,F_4,E_6,E_7`$, and $`E_8`$. Since nature is so exceptional, why not describe it using an exceptional Lie group?

Of the five exceptional Lie groups, four ($`G_2`$, $`F_4`$, $`E_7`$, and $`E_8`$) only have real or pseudoreal representations. A four-dimensional GUT model based on such a group will not give the observed chiral structure of weak interactions. The one exceptional group that does have complex or chiral representations is $`E_6`$, and this one works beautifully. The grand unified theory based on $`E_6`$ is not clearly superior to the $`SO(10)`$ model, but it does capture the successes of the $`SO(10)`$ model “exceptionally.”

But why $`E_6`$? The exceptional groups fit into a chain of embeddings

$$G_2F_4E_6E_7E_8.$$

If nature likes exceptional groups, why stop half-way? Yet $`E_6`$ is the only exceptional group that works for four-dimensional GUT’s.

There is another interesting chain of group embeddings (popularized notably by Olive),

$$SU\left(5\right)SO\left(10\right)E_6E_7E_8.$$

Here, at each step, one adds another node at one end of the Dynkin diagram. It is notable that the most significant GUT models correspond to the first three gauge groups in this chain. But again, we have the same question that we asked in connection with eqn. (1): if nature likes this chain, why stop half-way? As before, in four-dimensional GUT’s, we can only go half-way down this chain because of the $`VA`$ structure of weak interactions.

The next developments that I will mention were experimental. More accurate measurements showed that $`\mathrm{sin}^2\theta _W`$ is close to the GUT value, but not close enough. And the proton lifetime turned out to be longer than predicted in the simplest GUT’s.

Both of these problems were neatly addressed (by Dimopoulos, Raby, and Wilczek) by repeating the Georgi-Quinn-Weinberg calculation in the presence of supersymmetry. Supersymmetry raises the GUT prediction for $`\mathrm{sin}^2\theta _W`$, which becomes very close to the modern measurement. It also raises the GUT scale, making the proton lifetime long enough to be consistent with experiment, and lowering the gap between $`M_{GUT}`$ and the Planck mass $`M_{Pl}`$.

For these and other reasons, since the early 1980’s, SUSY-GUT’s have been the attractive form of GUT’s. One reason that I have not yet mentioned, which is very important even though it only involves a partial success, has to do with the hierarchy problem. SUSY stabilizes the hierarchy $`M_W/M_X<10^{13}`$, canceling large radiative corrections to the Higgs boson mass. A Higgs boson light compared to the GUT scale is thus made technically natural, but is not yet explained, and one is also left without an understanding of doublet-triplet splitting (the fact that the color triplet partners of the ordinary Higgs doublets are so heavy compared to the doublets; otherwise there would be very rapid proton decay). Also, it was shown (by Alvarez-Gaumé, Polchinski, and Wise) that radiative symmetry breaking would lead to a natural mechanism of electroweak symmetry breaking if the top quark is heavy enough. The requisite top quark mass, however, seemed bizarre at the time.

To continue the story, we must consider developments involving extra dimensions and gravity. It was soon realized that although one cannot unify three generations in four dimensions, one can readily do so if one starts above four dimensions. For example, I constructed an $`SO(12)`$ model in six dimensions and an $`SO(16)`$ model in ten dimensions, in each case getting, after compactification to four dimensions, three generations of quarks and leptons from a single irreducible representation of the unified group in higher dimensions.

A much bigger change came in 1984; with Green-Schwarz anomaly cancellation and the construction (by Gross, Harvey, Martinec, and Rohm) of the heterotic string, it became feasible to combine GUT’s with string theory and thus to unify all the forces, including gravity. In this framework (and assuming at least a small range of energies in which field theory ideas can be applied), one has to unify all the observed forces plus three generations of quarks and leptons plus Higgs bosons in one SUSY-multiplet – because that is all there is.

One also has to start in ten dimensions with $`E_8\times E_8`$ or $`SO(32)`$ because those are the only ten-dimensional gauge groups that are allowed by the anomaly cancellation mechanism. Of these two choices, the one that works is $`E_8\times E_8`$. (This assertion again assumes at least a small range of validity of field theory and was shown in early work on compactification by Candelas, Horowitz, Strominger, and me.)

So in incorporating gravity and string theory into the picture, one is forced to continue the GUT chain of eqn. (1) to the end, and to unify the three generations of quarks and leptons. We recall that these two steps did not work in four-dimensional GUT’s. The model is constructed by starting with $`𝐑^4\times K`$, where $`𝐑^4`$ is four-dimensional Minkowski space and, to preserve four-dimensional supersymmetry, $`K`$ is a compact six-manifold of a special sort, a “Calabi-Yau manifold.” Then, to obey the equations of motion, one is forced to introduce vacuum expectation values or VEV’s for gauge fields on $`K`$, breaking $`E_8`$ to a subgroup. This step works without any contrivance and is forced on us by the equations of motion.

By making a very simple choice (taking the gauge fields VEV’s to lie in an $`SU(3)`$, $`SU(4)`$, or $`SU(5)`$ subgroup of $`E_8`$), one can ensure that the unbroken subgroup of $`E_8`$ is one of the usual GUT subgroups $`E_6`$, $`SO(10)`$, or $`SU(5)`$. And the massless fermions occur in just the right representations, with the number of generations being a topological invariant of $`K`$ (together with its associated gauge bundle). But we do not know how to determine just what the gauge field VEV’s should be, and because many choices are possible for $`K`$, we also cannot predict the number of fermion generations.

These models are not really four-dimensional GUT’s, since unification really only occurs in ten dimensions. The predictions for fermion quantum numbers, $`\mathrm{sin}^2\theta _W`$, and proton decay are similar to those of four-dimensional GUT’s. There are also some differences. I will not survey the differences systematically, but I will mention a few of them. A higher-dimensional mechanism to solve the doublet-triplet splitting problem was pointed out in the original Calabi-Yau paper and used in many subsequent constructions. (This mechanism has been reconsidered in the last few years in the phenomenological literature and will be discussed in tomorrow’s lectures.) Also, in string theory models in which unification is achieved only in higher dimensions, the usual quantization of electric charge is generally not obeyed. There are superheavy unconfined particles with fractional electric charge, coming from strings wrapping around uncontractible loops in $`K`$. Perhaps they could serve as a dark matter candidate – though modern experimental bounds, from MACRO and elsewhere, put a severe limit on this possibility. Conversely, the quantum of magnetic charge is in these models generally larger than the Dirac quantum.

The other thing that happened in passing to string theory was that the gap between the GUT scale and the Planck scale was lowered again. We recall that in the original Georgi-Quinn-Weinberg computation, the GUT scale was below the Planck scale by a factor of about $`10^4`$. If we place things on a log scale, as we really should, then this is roughly a 20% discrepancy, and including supersymmetry reduced the discrepancy by about a third. Incorporating string theory, which for weakly coupled strings means roughly that the string scale enters the formula rather than the Planck scale, turned out to cut roughly another third off of the original discrepancy. So, by the mid-1980’s, the logarithmic discrepancy between the GUT and fundamental scales was cut to about 6%. What I think is important here is that the two improvements in the original computation – by including supersymmetry and string theory – both had independent justifications. They were not introduced just for this reason. The remaining discrepancy of 6% is small but real. There are interesting ideas about reducing it further or eliminating it (for example using the strongly coupled heterotic string instead of the weakly coupled one), but I think they are not as well motivated as were the enrichment of GUT’s by supersymmetry and string theory.

What has happened since? One important development is certainly that oscillations in solar and atmospheric neutrinos have pointed to neutrino masses (or at least mass differences) in roughly the range that had been estimated from GUT’s twenty years earlier. Meanwhile, astronomy has given other clues. The acceleration of cosmic expansion – whether ultimately interpreted in terms of a cosmological constant or a more elaborate mechanism perhaps involving a light scalar field – should eventually give an important clue about SUSY GUT’s. The small anisotropies that have been observed in the cosmic microwave radiation have for their simplest interpretation an early inflationary period of the universe at a scale near the SUSY GUT scale. So together with the neutrinos, this may represent two observations of phenomena involving an energy scale close to that of GUT’s, although alternative interpretations are possible in each case.

Another development is that the top quark mass turned out to be quite large – as assumed in supersymmetric models of electroweak symmetry breaking that were formulated in the early 1980’s. And, of course, increasingly precise tests of the standard model have been quite consistent with the SUSY-based approach to the hierarchy problem, while generally adding to the challenges faced by other approaches to that problem.

So in short, the GUT-based approach to physics has been attractive since it was first put forward close to thirty years ago; it has been enriched by new ideas, notably supersymmetry and strings; and there are real hints that it is on the right track, notably from $`\mathrm{sin}^2\theta _W`$ and neutrino masses.

If this approach is right, what may we find at accelerators? The Higgs boson really should be in reach, since grand unification does not work without an elementary Higgs boson, which moreover should not be too heavy or the standard model will break down at energies far below the GUT scale (spoiling the SUSY GUT prediction of $`\mathrm{sin}^2\theta _W`$). Just on this grounds, the Higgs boson probably should not weigh more than about 200 GeV, as is suggested in any case by the precision electroweak fits. If we take into account supersymmetry, which is also needed for grand unification, and assume the minimal supersymmetric spectrum, then the Higgs mass should be below about 130 GeV, and it is tempting to hope that it may lie at the 115 GeV value hinted at by LEP. If below 130 GeV, the Higgs should appear at Fermilab if the projected luminosity is reached, and in any case it should be seen at the LHC.

We can hope for much more beyond the Higgs. Supersymmetric particles should be in reach of the LHC – and maybe of Fermilab – since the supersymmetric approach to the hierarchy problem does not make sense if they are too heavy.

But what would the superworld really look like? Personally, I do not think that we have a convincing picture of this, in any detail, partly because of the problems in conservation of flavor, CP, and even baryon number. So even though we have a hunch that supersymmetry will be found, the details, if it is indeed found, will be surprising, at least to me.

The mystery about what the superworld might actually look like is part of what makes the search so exciting. Moreover, exploration of the superworld will be a long project, because there will be so many new particles and new interactions to unravel. It will require high precision, which hopefully will come from electron colliders such as TESLA or the LC, as well as requiring the extreme energy that can be obtained in proton colliders like Fermilab and the LHC. And interpreting the clues that will come from those machines will certainly require the best efforts of theorists.