text:
# Tycho Brahe’s supernova: light from centuries past

## 1 Introduction

SN 1572 was well observed in Europe (as well as in the Far–East) for almost two years. It added a new aspect to the debate at the time over the Aristotelian cosmological views, as it forced to reconsider the immutability of the heavens and the solid nature of the celestial spheres: the “star” gained brightness and lost it during a period of two years, but it showed no detectable parallax. According to prevalent views about the heavens, mutability would only happen in the sublunar region. This was even the place where comets were assumed to originate. The appearence of SN 1572 challenged the order of the celestial spheres. The observers who followed it up111The observers who mostly contributed to measure the position and luminosity, Tycho Brahe, Thomas Digges, Thaddeus Hagecius, Michael Mästlin, Jerónimo Muñoz, Caspar Peucer & Johannes Prätorius held very different views on the meaning of SN 1572. A comparison of their measurements and an account of their views is given elsewhere (Ruiz–Lapuente 2003). took sides with respect to established Aristotelian views. Today we can still use their observations to see whether that supernova would be of use for cosmology in a different way, as a distance indicator if seen by observers billions of years away from us. It is indeed possible to have a clear idea of where SN 1572 stands among its class.

After several centuries of questioning the nature of SN 1572, the identification as a Type I supernova came through the revision of the light curve done by Baade (1945) and based on data taken or quoted by Tycho Brahe (1603a). Before that time there were still speculations on whether it was a variable star of some kind, a nova (Morgan 1945) and it was even still considered the suggestion of its cometary nature (Lynn 1883). The cometary idea expressed in 1573 by Jerónimo Muñoz was based on the fact that the event is aligned with the Milky Way, so its decay in luminosity could be explained if it was a comet born among the stars that would first approach and then move away from us just along the line of sight. Tycho’s stella nova lies indeed only 49–98 pc above the Galactic plane.

In modern times, its comparison with other SNe allowed its classification as a Type I supernova (Baade 1945). Later on, this class was shown to contain events of very different nature: those identified with the explosion of white dwarfs (Type Ia) and those corresponding to the collapse of massive stars whose envelope had lost its hydrogen content in the interaction with a binary companion (Type Ib or c). SN 1572 was of the Type Ia class as discussed by de Vaucouleurs (1985) and by van den Bergh (1993). No further doubts that it is a Type Ia SN can now be held in view of the X–ray spectrum, which clearly differentiates SNIa events from SNe coming from the collapse of massive stars (Hughes et al. 1995). A first comparison of its luminosity at maximum with the bulk of SNeIa would have led to think that it was fainter than normal SNeIa. van den Bergh (1993) considers whether it could be a peculiar, subluminous SNIa, like SN 1991bg. However, SN 1572 was heavily obscured. It was reddened by E(B–V)=0.6$`\pm `$ 0.04 as it corresponds to the reddening of the stars near its position (Ruiz–Lapuente et al. 2003a,b). After taking into account the extinction undergone across the Milky Way, as well as the measurement of its decline rate, it is concluded that SN1572 was not a 91bg–like event. Neither was it an overluminous SNIa like SN 1991T or similars, but rather an event in the middle of the SNIa class.

## 2 New evidence on visual, color evolution and late decline

### 2.1  Visual light curve till 60 days

The date of maximum of SN 1572 is, from the available accounts, the most uncertain record. The first observation was done on Nov 6 though there is no record of the magnitude at that date (see for a full discussion Stephenson & Green 2002). On Nov 2 it was not yet noticeable as reported by Muñoz (1573). In Nov 11, date of the first observation done by Tycho Brahe (1603a) and the first reported observation by Muñoz (1573), it seems to have achieved a magnitude between those of Jupiter and Venus. The Nov 5 detection given in the report of the Progymnasmata by Tycho Brahe (1603b) is an erroneous date, as we have checked in the original. On Nov 16 and 17, when Caspar Peucer and Johannes Prätorius saw it, it seems to have a magnitude close to Venus (Table 1 and notes). Later on, in Jan 7, it was already fainter than Jupiter (Muñoz 1573). The reconstruction of the observations suggests an apparent visual magnitude at maximum of –4.0$`\pm `$ 0.3 (Baade 1945; de Vaucouleurs 1985). The overall light curve fitting done in the present work by trying different decline rates (see below) indicates that the maximum was broad and that it should have taken place around Nov 21. On Nov 11, and even on Nov 16, when Peucer and Prätorius observed it, it would be still on the rise. This is consistent with the mean risetime to maximum of SNeIa of 17.8$`{}_{1.0}{}^{}{}_{}{}^{+1.9}`$ days (Goldhaber et al. 2001). While most authors have taken arbitrarily Nov 15 as the date of visual maximum, it seems obvious that the maximum could not have happened before Nov 20.

The SN 1572 data are compared in this section to templates using the stretch factor s as characterization of the rate of decline (Perlmutter et al. 1999; Goldhaber et al 2001; Nobili et al. 2003). The stretch factor s method is used by the Supernova Cosmology Project to quantify the decline rate of the supernova from data extending up to 60 days after maximum. Even in absence of a measurement of the brightness at maximum, the method allows to produce a fine description of the supernova within the family of decline rates. The magnitudes and limits on magnitudes for SN 1572, when compared to templates, lead us to conclude that this supernova was not fast. The best agreement found for SN 1572 corresponds to s $``$ 0.9. In this regard, SN 1572 is very similar to SN 1996X which had a s = 0.889. A comparison can be seen in Figure 1, where the data of SN 1572 are shown overplotted to a s = 0.9 SNIa template in V. We also show for comparison the visual light curves of SN 1996X (normal SNIa of s $`=`$ 0.889) and SN 1991bg (a fast subluminous event of s $`=`$ 0.62). Therefore, these first 60 days show that there is no similarity with the fast subluminous SN 1991bg. A broad light curve as the one of SN 1991T can also be excluded from the constraints about its decline (see Figure 1).

### 2.2 Late–time decline

An additional proof that SN 1572 is a standard SNIa comes from its slow late decline between 100 and 450 days after maximum, similar to that of the bulk of SNeIa. SN 1572 had a late decline of 1.4 mag in 100 days during that epoch, in concordance with normal SNeIa and far from the fast decline of SN 1991bg (Filippenko et al. 1992; Leibundgut et al. 1993; Ruiz–Lapuente et al. 1993; Turatto et al. 1996). SN 1990N and other normal SNeIa had decline rates in such period of 1.38–1.5 mag in 100 days. In consistency with the findings near maximum, late–time and overall similarity to the light curves of normal SNIa is found. In order to test the adequacy of the fit using the whole light curve, we construct new templates for a normal event of s = 0.9 (well exemplified by SN 1996X), a subluminous 91bg–type event and an overluminous 91T–type event matching the available early and late–time observations. Before 60 days those templates are equal to the Hamuy et al. 1996b templates and the SCP templates for the corresponding stretches of those SNeIa. At later times they follow the available late–time photometry (Schmidt et al. 1993; Salvo et al. 2001) to their latest available epoch. The fit of the SN 1572 data to the SNIa template corresponding to s = 0.9 has a $`\chi ^2`$ of 14.44 for 10 degrees of freedom, which is acceptable. In contrast, the fit to the template of a fast declining, underluminous SNIa like SN 1991bg (s = 0.62) has an exceedingly high $`\chi ^2`$ of 53.55 for 4 degrees of freedom (the two premaximum points as well as the last 4 points have been excluded because of the lack of data for those subluminous SNeIa at such stage). In the other extreme, the fit to the template of a slow-declining, overluminous SNIa like SN 1991T (s = 1.2) is almost as bad, with a $`\chi ^2`$ of 82.99 for 10 degrees of freedom.

In Figure 2, we show the overall visual light curve till almost 500 days compared with that of SN 1996X and the templates for SN 1991bg and SN 1991T. The last upper limit given by Tycho Brahe indicates that the supernova had not yet leveled off at 480 days past maximum due to a light echo produced in an intervening dust cloud. Not all the SNeIa in dusty regions have shown light echoes associated with the scattering of the supernova light in neighboring dust clouds, though. The production of light echos depends on the spatial distribution of dust relative to both observer and supernova. In the case of SN 1986G, a very heavily reddened supernova, such light echo was not observed (Schmidt et al. 1993). However, for SN 1572 one can not exclude a leveling off occuring later, at a fainter magnitude level. SN 1991T and SN 1998bu have indeed shown a slowing down of the V magnitude decline at around 400 days as result of the light echo (Schmidt et al. 1993, Cappellaro et al. 2001), and the leveling off occurs at 500 days about 10 magnitudes below the maximum brigthness. For SN 1572 an equivalent step in the level of brightness would correspond to V$`=`$6, which is about the naked eye limit. SN 1998bu is a SNIa very similar to SN 1996X and SN 1572 but with a reddening of E(B–V)$`=`$ 0.32 $`\pm `$ 0.04 (Hernandez et al. 2000; Cappellaro et al. 2001). In fact, we can think of Tycho Brahe’ s supernova as SN 1998bu with twice the reddening.

Moreover, the late–time decline of SN 1572 informs us of deposition of the energy from the decay 56Co $``$ 56Fe in the ejecta of this supernova. Its slow decline would fit the declines predicted for SNeIa coming from the explosion of Chandrasekhar–mass WDs. The late–time decline in V indicates a significant deposition of the energy of the positrons which arise in the decay 56Co $``$ 56Fe. After 200 days since explosion, positrons become the dominant source of luminosity. Departures from full trapping of 56Co decay are seen in all SNeIa and they are expected. They arise from lack of confinement of positrons to the ejecta of the supernova and/or incomplete thermalization of their energies. Various effects can lead to the diversity in the late–time SN bolometric light curves: nucleosynthesis, kinematic differences, ejected mass, degree of mixing of the ejecta, and configuration and intensity of the magnetic field. In particular, as suggested in Ruiz–Lapuente & Spruit (1998), those late–time declines can serve to investigate the magnetic field configuration of the SN and trace it back to the original configuration of the WD right before explosion. It is found that departures from the full trapping of 56Co decay of the order of 10–15$`\%`$ at about 400 days can be explained by the distribution of radioactive material. Larger departures, of 30–40$`\%`$ or larger, have to be interpreted in terms of lack of magnetic field confinement of positrons, or even enhancement of the escape favored by a radially combed configuration of the magnetic field. In the case of SN 1991bg the very fast drop of the late light curve suggests not only a total ejected mass well below the Chandrasekhar mass, but an absence of a significant tangled magnetic field. An exploration of the late–time decline predicted for different explosion models compared with a sizeable sample of bolometric light curves was done by Milne, The and Leising (1999). For SN 1572, the availability of data is limited to the visual band. Though, analogies can be drawn with similar SNeIa for which the bolometric data are available.

The late decline is, as seen above, a major proof which confirms that this historical supernova is not part of the estimated 16$`\%`$ of intrinsically subluminous SNeIa (Li et al. 2001), likely arising from a different class of explosion ejecting a smaller amount of mass (Ruiz–Lapuente et al. 1993). If it were part of the class of peculiar subluminous SNeIa, caution should be given when used for cosmology, because it is not clear whether those subluminous events are in the linear relation brightness–rate of decline.

### 2.3 Reddening and color

The most direct estimate of the reddening comes from our measurement of the reddening and extinction of the stars around the geometrical center of SN 1572 at the distance of the supernova (Ruiz–Lapuente et al. 2003a,b). The stars which are compatible with a distance in the range of 3$`\pm `$1 kpc (see de Vaucouleurs 1985 for a review of the distances towards SN 1572) have reddenings which average to E(B–V)=0.6 $`\pm `$ 0.04. Those values were derived in the programme of search of the companion star of Galactic SNeIa (Ruiz–Lapuente et al. 2003a,b). Star candidates within a radius compatible with the degree of precision on the center and distance gained during the 431 years since explosion are modeled to determine Teff, log g, metallicity, distance and E(B–V). The radial velocities are measured. The individual reddening values found in that programme for the stars within 35$`\%`$ of the radius of the remnant go from E(B–V) = 0.50 to E(B–V)=0.8, the values increasing with distance.

The above mean reddening to SN 1572 corresponds to an extinction $`A_V=`$ 1.86 $`\pm `$ 0.12 mag if we take $`R_V=`$ 3.1 (Sneden et al. 1978; Riecke & Lebofsky 1985). The Galactic extinction data based on COBE/DIRBE observations (Drimmel & Spergel 2001) give a value of $`A_V=1.77mag`$ in that direction (the maximum Galactic extinction being $`A_V=1.90mag`$ there). Thus the extinction measured at the distance of SN 1572 is in agreement with the COBE results. Correcting the apparent brightness from the dimming by dust, we find again that Tycho’s SN was not subluminous but a “normal” SN Ia.

We can correct as well the color of the supernova from the records of the epoch and derive the intrinsic color evolution. SNeIa are found to show a similar color evolution two months after maximum, with a well–established law of low intrinsic dispersion valid from 30 to 90 days, explored in Phillips et al. 1999 (Ph99) based on the work of Lira (1995):

$$(BV)_0=0.725(\pm 0.05)0.0118(t_V60)$$
$`(1)`$
where $`t_V`$ is the time since visual maximum.

The color of SN 1572 two months after discovery was reported to be similar to Mars and Aldebaran (thus B–V in the range 1.36–1.54). Previous and subsequent color estimates are also shown in Table 2. After correcting the observed color by our measurement of reddening $`E(BV)=`$ 0.6 $`\pm `$ 0.04, the intrinsic color at 55 $`\pm `$10 days is $`(BV)_0=`$ 0.76 $`\pm `$ 0.24, very much in agreement with the expected $`(BV)_0=`$ 0.78 $`\pm `$0.15 for that epoch. Such color fits very well in the Ph99 law for the time evolution of the intrinsic color of the SNeIa. In agreement with that law, Nobili et al. (2002) have shown that the color evolution of the bulk of SNeIa has a low dispersion of 0.1 mag in the tail. Before maximum the corrected color of SN 1572 would be $`(BV)_0=`$ 0.22 $`\pm `$ 0.29, thus consistent with the fact that normal SNeIa have $`(BV)_0`$ $``$ 0. SN 1991bg and other subluminous SNeIa clearly deviate from the standard color evolution at early epochs, being intrinsically red at maximum $`(BV)_0=0.6`$ (Leibundgut et al. 1993; Ph99). Moreover, the color evolution derived from the B and V templates for a s $`=`$ 0.9 SNIa gives a good agreement with the supernova. In Figure 3 we show the color light curve of SN 1572, corrected by its E(B–V)$`=`$ 0.6, compared with that of “normal” SNeIa. Both the peculiar SN 1991bg and the normal SN 1996X are also shown in their reddening–corrected colors.

The data displayed in this figure come from different observers (see Table 2 and notes). Color estimates by eye are difficulted by the lack of an established system which encompases the whole range of colors. The conversion of eye–estimated colors to our filter scale is addressed in Pskovskii (1977) and Schaefer (1996). To add to the inherent eye uncertainty in estimating color, it has to be noted that the first comparisons done by the observers of the time using planets like Jupiter, Venus or Saturn as references, were difficulted by the fact that the planets were often very far from the position of SN 1572 (away in the sky by 590). De Vaucouleurs (1985) discussed the effect of that on the early visual estimates (later estimates since they are done through comparison with particular stars are more reliable).

In the nebular phase, at 175 days according to Tycho Brahe’s records (Table 2) the supernova went back to a white color, which is consistent, after correcting for extinction, with what is found in normal SNeIa as well. Uncertainties in the color estimates have been assigned and new records included.

### 2.4 The luminosity of Tycho’s supernova

Given the apparent visual magnitude, the luminosity of Tycho Brahe’s supernova would have corresponded to a visual absolute magnitude of $`M_V=17.725log(d/3.5kpc)A_V`$. If corrected for $`A_V`$ = 1.86 $`\pm `$0.12, this means $`M_V=19.585log(d/3.5kpc)\pm 0.42`$.

De Vaucouleurs (1985) gives as the most likely estimate to the distance to SN 1572 3.2 $`\pm `$ 0.3 kpc. Taking 3.2 kpc for the distance, the luminosity would be –19.38 $`\pm `$ 0.42. That compares well with the value $`M_V=19.12\pm 0.26`$, the mean magnitude from the Calán/Tololo sample (Hamuy et al. 1996a).

To get a more refined distance to SN 1572, one would need to follow a similar procedure as in Winkler, Gupta & Knox (2003) for SN 1006. These authors obtained a distance to SN 1006 of 2.18 $`\pm `$ 0.08 kpc using the SNR expanding parallax method (Kirshner, Winkler & Chevalier 1987) with precise measurements of the proper motion of the filaments in SN 1006 and narrower estimates of the shock velocity of the expanding remnant.

The proper motion of the filaments in SN 1572 was measured by Kamper & van den Bergh (1978), and it is about 0.22 $`\pm `$ 0.02 $`^{^{\prime \prime }}`$ yr-1. There is, however, a larger uncertainty in the estimate of the shock velocity, which is obtained through modeling of the broad H$`\alpha `$ emission of the material swept up by the SN. The estimate is model dependent and improves with a better knowledge of the composition of the fully ionized plasma behind the shock. Kirshner, Winkler & Chevalier (1987) had suggested from the observations available at that time a distance to SN 1006 between 1.4 and 2.1 kpc and a distance to SN 1572 of 2.0–2.8 kpc. The estimates of the shock velocity, though, could soon be revised from the wealth of spectral observations on SN 1572’s filaments taken during the last years.

A rough estimate of the distance to Tycho’s SN can be made by comparison of the angular diameters of the two historical SNeIa. The angular diameter of G327.6+14.6 (the remnant of SN 1006) is $`\theta (1006)=30.60\pm 0.06`$ arcmin whereas that of Tycho is $`\theta (Tycho)=8.65\pm 0.05`$ arcmin. We can place $`\theta (Tycho)/\theta (1006)=d(1006)/d(Tycho)t(Tycho)/t(1006)\eta `$, where $`t(Tycho)`$ and $`t(1006)`$ are the respective ages of the remnants and $`\eta =v(Tycho)/v(1006)`$ is the ratio of the average expansion velocities of the two SN ejecta. Such approach gives, if we use the distance obtained by Winkler, Gupta & Long (2003) to SN 1006, a distance to SN 1572 of $`d(Tycho)3.3`$ kpc, well in the middle of other estimates.

All of the above shows that it is at present hard to escape from the uncertain range in the distance to SN 1572, i.e. from 2 to 4 kpc with the analyses done so far. Such limitation does not prevent to infer the right place of SN 1572 within the luminosity sample: the conclusion on the brightness–rate of decline for SN 1572 is based on the gathered empirical evidence and the well–established use of the peak luminosity–decline correlation.

Incidentally, if a good distance would be available for SN 1572, one could obtain the distance to SN 1996X, which is not well determined with the present available methods. The Dn$`\sigma `$ relation gives a distance modulus $`\mu `$ = 31.32 $`\pm `$ 0.4 (Faber et al. 1989), whereas the SNeIa calibration using H0 = 65 km s-1 Mpc-1 would place it a significantly larger distance (Salvo et al. 2001). The very good agreement between the late time light curves of SN 1572 and SN 1996X allows us to infer a $`\mathrm{\Delta }\mu `$(SN1572–SN1996X) $`=`$ 18.9 $`\pm `$ 0.12 between both SNeIa.

## 3 Conclusions

We have found that SN 1572 was a supernova very close to the template with a stretch factor s $``$ 0.9 (see Table 3 for a summary). The light curve grows in precision towards the late times, being highly uncertain around maximum brightness. An overall agreement between early, late decline and color with the expected evolution of normal SNeIa supports our conclusion.

Type Ia supernovae with stretch factors between 0.9 and 1.1 make the vast majority of the observed population. They are not only those most frequently found in nearby searches, but also the bulk of discoveries in cosmological searches at high–z, as can be seen in the sample of SNe at z $`>`$ 0.3 found by the Supernova Cosmology Project (Perlmutter et al. 1999). Among SNeIa of s $``$ 0.9 in nearby galaxies for which very late–time data are available, we have found a close resemblance to SN 1996X in rate of decline. SN 1572 likely has a slightly slower rate. However, whereas SN 1996X was not heavily reddened, the reddening in SN 1572 is E(B–V)=0.60 $`\pm `$ 0.05. While at present heavily reddened SNe Ia involve a larger difficulty than unreddened ones for their proper use in cosmology (Riess et al. 2000; Knop et al. 2003), a better use of multicolor light curves should allow to overcome the problem in their calibration for cosmology.

Tycho Brahe’s supernova belongs to the SNIa type that are used for the determination of the nature of dark energy. Its light, traveling already light years away from us, carries the information on our acceleration epoch to an even more dark–energy dominated future. If used for cosmology by observers in galaxies billions of light years distant, it would require, certainly, a careful handling of its intrinsic brightness.

Acknowledgements This paper amongst others to come is dedicated to my father. Past, present and future merged for him just a few months ago, as sadly as unexpectedly. His life was devoted to help the souls of many and to carry research in brain and mind studies. His inner light lives forever with us. I would like to express my gratitude to M. A. Granada, from the Department of History of Philosophy of the U. Barcelona, for helping to locate the XVI century historical sources.

NOTES TO TABLE 1

1 Jerónimo Muñoz writes that he was certain that on November 2, 1572 the “comet” was not in the sky, as he was teaching to know the stars to a numerous group of people that evening (Muñoz 1573).222 Jerónimo Muñoz (?-1592) was professor of Hebrew and Mathematics at the University of Valencia and later on of Astrology at the University of Salamanca. The records assembled here are from his Book of the New Comet (Libro del nuevo cometa) written as King Philip II consulted him on SN 1572. He determines the absence of parallax of the star and gives estimates of color and magnitudes, though his account is mostly addressed to criticize the Aristotelian cosmology. He wrote that the heavens and stars are not made of quintessence, but related to the elements, and advocated to reconsider the presocratic cosmological views of Democritus and Anaxagoras. Tycho Brahe (1603b) compiled parts of Muñoz’s observations and discussed those views in the third volume of the Progymnasmata. He left out, however, some useful information that we have added here. We corrected as well the premaximum quote.

2 Nov 11 is the first observation quoted by Muñoz in his book. The apparent magnitude appeared somewhat greater than that of Jupiter, which was 59 degrees distant from the comet, and almost equal to that of the Morning Star. Jupiter was at that time V=–2.61 and Venus V=–4.35 (Baade 1945). We assign –3.0$`\pm `$0.3.

3 On Nov 11 Tycho Brahe also saw it, and his preliminar report in De nova stella said that equaled Venus when this planet was at maximum brightness. He was, however, more conservative in the account he gives in Progymnasmata. We give –3.5 mag to the estimate of Tycho $`\pm `$0.3. In the following we must keep in mind that Tycho Brahe observes from a latitude of $``$ N550, as he stays in the area of Copenhague for two years since the discovery of SN 1572. From that location the supernova can be observed well above the horizon through the whole year. De Vaucouleurs notes that when comparing SN 1572 with Venus one has to account from the fact that the earlier is not affected by extinction while the second was significantly affected. This places V=–3.8 or fainter if one considers the apparent magnitude of Venus affected by extinction.

4 Tycho Brahe reports in the Progymnasmata on the observations done by Caspar Peucer and Johannes Prätorious on Nov 16. According to the first, it surpassed in brightness all stars and planets with the exception of Venus, for the second, it was brighter than Jupiter but easily fainter than Venus. Michael Mästlin who saw it before Nov 17 states that almost equalled Venus. It has been given to these estimates a value of V=–4.0 $`\pm `$0.3.

5 Tycho Brahe does not give the exact date of the observations. He quotes the month or an interval encompassing two months. It is reasonable to think that he refers to dates in the middle of the period. Therefore, we follow here Baade (1945) and assign the middle of the period as date of observation. However, we add $`\pm `$ 7 days of uncertainty in the case a month is quoted, and $`\pm `$ 14 days if the period encompasses two months.

6 This is an upper limit given by Muñoz. He writes, when it began to be visible, it looked larger than Jupiter, and now, on January 7, 1573, it already looks smaller than Jupiter; this could have happened because it had risen higher than where it was when it first appeared. This upper limit to the rate of decline as well as the data around it by Tycho Brahe constrain the light curve stretch factor below s $`=`$ 1.1.

7 This is a difficult estimate. The interval quoted encompasses 2.46 mag (see Baade 1945). It is uncertain the interpolation done, from Jupiter with V=–2.18 to the brightest stars of the first class which are +0.28. We assign to the Baade estimate of V=–1.4 an uncertainty of $`\pm `$ 0.4.

8 This a bit better estimate, as encompasses an smaller range in magnitudes $`\pm `$ 0.2 if one selects the brighter range of first magnitude stars (0.2–0.3).

9 It is a very ambiguous quote, one should give it $`\pm `$ 0.5. The stars of second magnitude had a range encompassing 1 mag, $`m_V=`$4.0 and $`\sigma _V=`$0.5.

This is the point that lies far from the expected curve. It has been checked whether a problem of differential extinction could have affected SN 1572 with respect to the second magnitude stars, but there were second magnitude stars at similar altitude as SN 1572. The offset of 0.6 mag can not be attributed to extinction problems, but to ambiguous description of the record.

10 Tycho Brahe’s estimates in magnitude are more precise now as he refers to particular stars not far from the supernova with a low scatter in magnitude. The mean of those stars is $`m_V=`$2.48 and $`\sigma _V=`$0.23. We assign $`\pm `$0.2 to the error in the mags of those stars.

11 It becomes a very precise estimate by comparing to the fourth mag stars in Cassiopeia, when knowing that it was brighter than the nearby eleventh star of Cassiopeia. An error of 0.2 mag is assigned.

12 This seems to be the most precise description as it says it is equal to $`\kappa `$ Cas.

13 Errors of 0.2 and 0.4 mag respectively seem reasonable for these last estimates.

14 It can be noted that the quote that it was not observable later than March 1574 (V $`>`$ 6.0) allows to exclude that the luminosity levelled off, due to a light echo, at $``$ 10 mag below the supernova maximum brightness.

NOTES TO TABLE 2

1Various reports (see Baade 1945) give a color similar to the colors of Venus and Jupiter.

2 Tycho Brahe (1603a) states that after maximum it became yellowish in color. By comparison to the description as whitish (like Jupiter) and red like Aldebaran, we assign B–V $`=`$1.0 $`\pm `$0.25.

3 Muñoz (1573) states that the color was between those of Saturn and Mars, closer to Mars. One could identify then the color as 1.2 $`\pm `$ 0.1. Unfortunately, there is ambiguity in the date he is refering to: he mentions that he first saw it on Dec 2, and checked with observers that it was already there on Nov 11. He uses the expression “the color then was”, so it is not clear whether he refers to Nov 11 or Dec 2. It seems that the description, because of its detail, is based on his own appreciation of color on Dec 2 (we place an errorbar in the dates towards Nov 11, and a 7 days error after Dec 2).

4Tycho Brahe notes that shortly after the beginning of the year it was red like Aldebaran (Brahe 1603a). We assign 0.25 mag of uncertainty to those estimates.

5 According to Prätorius it was like Mars two months after its appearance (around Jan 15).

6Mästlin refers that at the end of February it is going back to its whitish color. Tycho states that it had a red color along February and March, becoming whitish after May. However, this seems to be a worse account of the color than the one given by Mästlin. In April it is described as silverish by some observers in Spain, and as similar to Saturn (van den Bergh 1970).

7 The last estimate is given 0.82 $`\pm `$0.25 in account for Tycho Brahe on the color in May.

Normal SNe Ia of s = 0.9. Whereas the $`\mathrm{\Delta }m(V)_{20}`$ is not well defined for SN 1572 due to a lack of accurate observations in the week around the maximum, the rate of decline in terms of stretch can be well determined by fitting premaximum and postmaximum up to 60 days.


format: text

text:
# Occupation number and fluctuations in the finite-temperature Bose-Hubbard model

## Abstract

We study the occupation numbers and number fluctuations of ultra-cold atoms in deep optical lattices for finite temperatures within the Bose-Hubbard model. Simple analytical expressions for the mean occupation number and number fluctuations are obtained in the weak-hopping regime using an interpolation between results from different perturbation approaches in the Mott-insulator and superfluid phases. These analytical results are compared to exact one dimensional numerical calculations using a finite temperature variant of the Density-Matrix Renormalisation Group (DMRG) method and found to have a high degree of accuracy. We also find very good agreement in the crossover “thermal” region. With the present approach the magnitude of number fluctuations under realistic experimental conditions can be estimated and the properties of the finite temperature phase diagram can be studied.

The Bose-Hubbard (BH) model, first studied in the late 1980’s Fisher89 , has long been considered a rather academic testing ground for analytical as well as numerical approaches in quantum statistical mechanics. Its most interesting feature is the existence of a quantum phase transition Sachdev  at zero temperature between a superfluid and a Mott-insulating phase. With recent advances in the experimental techniques of atom trapping and optical lattices, the BH model has regained substantial practical importance. As first noticed in Jaksch98 , the BH Hamiltonian applies to bosonic atoms trapped in a deep lattice potential. Recent experiments with optical lattices Greiner02  have spectacularly confirmed the existence of the superfluid-insulator transition. Since the Mott-phase is characterized by a well defined occupation number of the potential wells, the phase transition has important potential applications in quantum information processing Bloch-QI-exp ; Briegel , or Heisenberg-limited matter-wave interferometry Kasevich , both of which require optical lattices with regular filling. Although techniques to eliminate lattice defects have been designed Cirac  an important issue for the implementation of these applications under experimental conditions are fluctuations in the particle number per site at finite $`T`$. In this work we analyze these fluctuations by studying the finite temperature Bose-Hubbard model.

The zero temperature BH model has been extensively investigated in the past by methods such as the Gutzwiller projection ansatz RokhsarKotliar91 , the strong-coupling expansion FreericksMonien94 ; FreericksMonien96 ; ElstnerMonien99A ; ElstnerMonien99B  and Quantum Monte-Carlo KrauthTrivedi91 ; BatrouniScalettarZimanyi90 ; Scalettar95 ; Batrouni02 ; Svistunov02 , as well as with various mean-field approaches (see, e.g., Fisher89 ; Amico98-00 ; Stoof01 ). A very powerful numerical technique in the case of one spatial dimension is the density matrix renormalization group (DMRG) DMRG , used to calculate first and second order (i.e., amplitude and number-number) correlations KuenerWhiteMonien00 . Less attention has been paid to the nonzero temperature properties of the BH system. We note here that the quantum phase transition exists in the strict sense only at $`T=0`$. For finite $`T`$ the compressibility, i.e. the change of the particle number per site with the chemical potential, is always nonzero and thus only an approximate Mott phase exists. Furthermore the transition between superfluid and insulating behaviour passes through an intermediate thermal region. Recently this thermal crossover has been studied within a slave-boson formalism Stoof03 . In this work we study the finite-$`T`$ properties using a perturbative analytic approach in the strong-coupling limit, deriving simple analytical expressions for the occupation number and number fluctuations which cover both quantum and thermal effects. To verify these results numerically in one spatial dimension, we employ a finite-$`T`$ version of the DMRG approach DMRG . We observe good agreement between perturbative and DMRG results, including the description of the thermal crossover between the Mott and insulating phases.

We consider a $`d`$-dimensional cubic lattice of $`N`$ nonlinear oscillators characterised by the following Hamiltonian,

$`\begin{array}{c}\hfill \widehat{}\mu \widehat{𝒩}=\widehat{}_{\text{hop}}+\widehat{}_0,\end{array}`$ (1)
where $`\mu `$ is the chemical potential, $`\widehat{𝒩}=_𝒌\widehat{n}_𝒌`$ is the operator of the total number of atoms, and $`\widehat{n}_𝒌=\widehat{a}_𝒌^{}\widehat{a}_𝒌`$, with $`\widehat{a}_𝒌^{}`$, $`\widehat{a}_𝒌`$ being the usual bosonic creation/annihilation pairs. The lattice sites are numbered with a $`d`$-dimensional index $`𝒌=\{k_1,\mathrm{},k_d\}`$, while $`_𝒌`$ means summation over all lattice sites. The first term,

$$\widehat{}_{\text{hop}}=J\underset{𝒌𝒍}{}\left(\widehat{a}_𝒌^{}\widehat{a}_𝒍+\widehat{a}_𝒍^{}\widehat{a}_𝒌\right),$$
(2)
on the RHS of (1) is the hopping Hamiltonian, $`J0`$; the sum extends over nearest neighbours and we assume cyclic boundary conditions. The second term on the RHS of (1), $`\widehat{}_0=_𝒌H_{\text{nl}}\left(\widehat{n}_𝒌\right)`$, is a sum of the nonlinear site Hamiltonians

$$\begin{array}{c}H_{\text{nl}}\left(\widehat{n}_𝒌\right)=\frac{U}{2}\widehat{n}_𝒌\left(\widehat{n}_𝒌1\right)\mu \widehat{n}_𝒌\hfill \\ \hfill =\frac{U}{2}\left(\widehat{n}_𝒌\overline{n}\right)^2+\text{const},\end{array}$$
(3)
where $`\overline{n}=\mu /U+1/2`$ is the filling. The quantum averaging is defined in the standard manner as $`\mathrm{}=Z^1\text{Tr}\left[\text{e}^{\beta \left(\widehat{}\mu \widehat{𝒩}\right)}(\mathrm{})\right]`$, where $`Z=\text{Tr}\text{e}^{\beta \left(\widehat{}\mu \widehat{𝒩}\right)}`$ is the statistical sum, and $`\beta =1/T`$ is the inverse temperature. Note that oscillator units with $`\mathrm{}=k_\text{B}=1`$ are used throughout this paper.

While elaborate methods are needed if we wish to exactly pinpoint lobe boundaries or calculate long-range correlations, local quantities like the average number of atoms and the number fluctuations may be obtained at lesser cost. Let us firstly consider the Mott insulator phase. To zero order in hopping, the ground state of the Mott insulator is $`|0_{J=0}=_𝒌|n_0_𝒌,`$ where $`n_0=\text{round}(\overline{n})`$, cf. Eq. (3). A first order perturbative correction to this state will contain all states which are found from $`|0_{J=0}`$ by moving one atom to a neighboring site. All such “single-hop” states are eigenstates of $`_0`$ with the same relative energy $`U`$, resulting in an especially simple expression for the ground state to first order in $`J`$ Burnett03 . Up to second order we find

$`\begin{array}{c}\hfill |0=\alpha \left(𝟙{\displaystyle \frac{\mathrm{𝟙}}{𝕌}}_{\text{hop}}\right)|0_{J=0}+𝒪\left(_{\text{hop}}^2\right),\end{array}`$ (4)
where $`\alpha `$ is a normalisation factor. By direct calculation we then find

$$\begin{array}{c}0\left|\widehat{n}_𝒌\right|0=n_0,\\ \delta n^2=0\left|\widehat{n}_𝒌^2\right|00\left|\widehat{n}_𝒌\right|0^2=\frac{2zJ^2n_0(n_0+1)}{U^2},\end{array}$$
(7)
where $`z=2d`$ is the number of nearest neighbours in the lattice. These formulae are expected to be a good approximation deep in the Mott insulator phase where a large energy gap exists and the thermal occupation of higher levels can be disregarded.

It is worth noting that, if staying within the perturbation approach, these expressions cannot be extended to the thermodynamic limit $`N\mathrm{}`$. To demonstrate this, consider the formula thus obtained for the normalisation constant, $`\alpha `$:

$$\begin{array}{c}\alpha ^2\left(1+\frac{1}{U^2}0\left|_{\text{hop}}^2\right|0\right)=1,\\ 0\left|_{\text{hop}}^2\right|0=zJ^2Nn_0(n_0+1).\end{array}$$
(10)
These results hold only if $`zJ^2Nn_0(n_0+1)U^2`$, therefore making the thermodynamic limit impossible. (This problem went unnoticed in Ref. Burnett03 .) On the other hand, local (on-site) quantities should become independent of the size of the system when it becomes large enough (thermodynamic limit). In fact, our DMRG simulations show that Eq. (7) equally holds if $`zJ^2Nn_0(n_0+1)U^2`$.

In the superfluid phase, one has to account for states with the total number of quanta different from $`Nn_0`$. In the limit of zero hopping an arbitrary eigenstate is a product state, $`|\stackrel{}{n}=_𝒌|n_𝒌_𝒌`$. With $`\overline{n}`$ close to a half-integer, it suffices to keep only two states for each site, $`|n_0`$ and $`|n_0+1`$, effectively turning bosons into fermions. We shall call the subspace of states with $`n_𝒌=n_0,n_0+1`$, the f-subspace. In this subspace we introduce standard fermionic creation and annihilation operators, $`\widehat{c}_𝒌\widehat{c}_𝒍^{}+\widehat{c}_𝒍^{}\widehat{c}_𝒌=\delta _{𝒌𝒍}`$. Within the f-subspace, ($`l,m0`$)

$`\begin{array}{cc}\hfill \left(\widehat{a}_𝒌^{}\right)^{l+m}\widehat{a}_𝒌^l& =\left(c_𝒌^{}\sqrt{n_0+1}\right)^m\left(n_0+\widehat{c}_𝒌^{}\widehat{c}_𝒌\right)^l,\hfill \\ \hfill \left(\widehat{a}_𝒌^{}\right)^l\widehat{a}_𝒌^{l+m}& =\left(n_0+\widehat{c}_𝒌^{}\widehat{c}_𝒌\right)^l\left(c_𝒌\sqrt{n_0+1}\right)^m.\hfill \end{array}`$ (11)
Using these relations one can express any bosonic operator in fermionic terms. In particular,

$`\begin{array}{c}\hfill H_{\text{nl}}={\displaystyle \frac{U}{2}}\left(\widehat{c}_𝒌^{}\widehat{c}_𝒌+n_0\overline{n}\right)^2=\mathrm{\Delta }E\widehat{c}_𝒌^{}\widehat{c}_𝒌+\text{const},\end{array}`$ (12)
where $`\mathrm{\Delta }E=U\left(1/2+n_0\overline{n}\right)`$. Projection onto the f-subspace thus makes the Hamiltonian linear so it can be directly diagonalised: (omitting an additive constant)

$$\widehat{}\mu \widehat{𝒩}=\underset{𝒍}{}\left(E_{\phi _{l_1}}+\mathrm{}+E_{\phi _{l_d}}\right)\widehat{f}_𝒍^{}\widehat{f}_𝒍,$$
(13)
where $`\phi _l=2\pi (l1)/N`$, $`E_\phi =\mathrm{\Delta }E2J(n_0+1)\mathrm{cos}\phi `$, and $`\widehat{f}_𝒍`$ are related to $`\widehat{c}_𝒌`$ by a $`d`$-dimensional discrete Fourier transform. The average on-site particle number and number fluctuations are then expressed by the average number of “fermions” per site, $`\mathrm{\Delta }n=\widehat{c}_𝒌^{}\widehat{c}_𝒌`$, as

$`\begin{array}{cc}\hfill \widehat{n}_𝒌& =n_0+\mathrm{\Delta }n,\hfill \\ \hfill \widehat{n}_𝒌^2& =\left(n_0+\widehat{c}_𝒌^{}\widehat{c}_𝒌\right)^2=n_0^2+(2n_0+1)\mathrm{\Delta }n,\hfill \\ \hfill \delta n^2& =\widehat{n}_𝒌^2\widehat{n}_𝒌^2=\mathrm{\Delta }n(1\mathrm{\Delta }n).\hfill \end{array}`$ (14)
For $`\mathrm{\Delta }n`$ we readily find: ($`\beta =1/T`$)

$`\begin{array}{c}\hfill \mathrm{\Delta }n=\left({\displaystyle \underset{s=1}{\overset{d}{}}}{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\phi _s}{2\pi }}\right){\displaystyle \frac{1}{1+\mathrm{exp}\beta _{s=1}^dE_{\phi _s}}}.\end{array}`$ (15)
This formula holds under the thermodynamic limit $`N\mathrm{}`$.

The relations (7) and (14), (15) are found under mutually exclusive conditions. Equation (7) follows if we assume that only the ground state matters, and retain the first non-vanishing correction to its wave function; with only one level being important temperature is of no concern. Conversely, equations (14), (15) ignore corrections to the eigenstate wavefunctions yet account for energy shifts; with many levels being accounted for, thermal properties are retained. The question we now address is whether a relation can be derived which covers both insulator and superfluid regions, plus, even more importantly, the crossover region (which is naturally termed the thermal region)? It would be natural to include perturbative corrections to all “fermionic” states, but this results in some intractable algebra. It turns out, however, that for $`JU`$, simple interpolating formulae betwen the two above results may be found which apply not only to the insulator and superfluid regions, but also to the thermal crossover region between them. Namely,

$`\begin{array}{cc}\hfill \widehat{n}_𝒌& =n_0+\mathrm{\Delta }n,\hfill \\ \hfill \delta n^2& =\mathrm{\Delta }n(1\mathrm{\Delta }n)+{\displaystyle \frac{2zJ^2(n_0+\mathrm{\Delta }n)(n_0+\mathrm{\Delta }n+1)}{U^2}},\hfill \end{array}`$ (16)
where $`\mathrm{\Delta }n`$ is given by Eq. (15). The idea behind (16) is relatively simple: for small $`J`$, the quantum contribution (7) is negligible in the thermal and (even more so) the superfluid regions, while in the insulator region, $`\mathrm{\Delta }n0`$ (or 1), so that (16) coincide with (7).

To verify these results in one spatial dimension, we use a finite-temperature version of the DMRG approach DMRG . At the core of our approach is the following block-doubling algorithm: assuming that we know the optimised states for a block of length $`L`$, we can then find eigenstates of a ring of three such blocks, (recall that we use periodic boundary conditions) and calculate the thermal $`\rho `$-matrix. Tracing over the states of one block yields the $`\rho `$-matrix of a double-sized block. Diagonalising this $`\rho `$-matrix and taking a certain number of eigenstates corresponding to the largest eigenvalues results in the optimised basis in a block of length $`2L`$. The loss of probability is assessed by summing the neglected eigenvalues; if this loss becomes unacceptable the iterations are stopped. To initiate the algorithm, we start from a ring of $`L_0+L^{}`$ oscillators, and then trace out the states of $`L^{}`$ oscillators, resulting in an optimised basis for a block of $`L_0`$ sites. We found this to yield much better results than starting from an open-ended block by simply selecting its lowest eigenstates. The full details of the algorithm will be discussed elsewhere.

A word of caution is necessary here. While borrowing much of the technical side of the DMRG method, our approach is conceptually different, the difference residing in the importance of temperature. Similar to other renormalisation-group methods, DMRG relies on the fact that the dimension of the physically relevant subspace of the full Hilbert space does not grow with the system size. For finite temperatures this can only be true until the thermal correlation length is reached. Consequently we only grow the block up to a size comparable to the thermal correlation length, and must simply stop as soon as the probability loss becomes unacceptable. We use periodic rather than open boundary conditions based on the same argument. For the temperatures considered here the maximum block length reached in this way is long enough such that finite-size effects are unimportant. For higher temperatures a combination of stochastic and DMRG techniques can be used, which we are developing and will be discussed in a later work.

The results of our calculations using both DMRG and perturbative approaches in 1-D are summarised in the figure, where we plot average on-site numbers and number fluctuations for three “cross-sections” along the $`\mu `$-axis, for $`J=0.01U,0.02U,0.05U`$. Each plot shows data for two temperatures, $`T=0.01U`$ and $`T=0.001U`$. The top and bottom rows of plots represent, respectively, the on-site population, $`\widehat{n}_k`$, and number fluctuations, $`\delta n^2=\widehat{n}_k^2\widehat{n}_k^2`$. To gain more insight into the thermal region, the middle row of plots shows the difference, $`\mathrm{\Delta }n`$, between the on-site population and the nearest integer value. At $`T=0.001U`$, the insulator region is clearly defined by $`\mathrm{\Delta }n`$ abruptly falling to zero. At $`T=0.01U`$, $`\mathrm{\Delta }n`$ is a smooth function vanishing as $`\mu `$ goes deeper into the insulator region. A similar effect is seen for the on-site number fluctuations. At $`T=0.001U`$, the phase transition points are well defined, whereas at $`T=0.01U`$ the boundaries of the insulator phase are “eroded” so that $`\delta n^2`$ is a visually smooth function of $`\mu `$.

Comparing the DMRG to perturbative results, we see that the latter provide a surprisingly good description of the quantities in question. For $`J=0.01U`$ (left column of plots), we find a very good agreement between the perturbative and DMRG results. This agreement is reasonable for $`J=0.02U`$, (center column) becoming only fair at $`J=0.05U`$ (right column). Note that, even in the latter case, the error is mostly in positioning the thermal region between the insulator and superfluid phases. Away from the thermal region, the results of the perturbative approach show good agreement with the DMRG results, both for the superfluid and for the insulator.

###### Acknowledgements.

This research was supported by the Deutsche Forschungsgemeinschaft under contract FL210/12 and the New Zealand Foundation for Research, Science and Technology (Grant No. UFRJ0001).


text:
# Human housekeeping genes are compact

## I Assignment of housekeeping genes

A recently published database provides microarray expression data for Affymetrix U95A chip, containing 12,600 probes, and hybridized to 101 different samples 9  from 47 different human tissues and cell lines. These samples are mainly from the normal human physiological state, and therefore this dataset provides a description of the normal mammalian transcriptome.

We calculated the distribution of the number of different tissues in which a gene is expressed. Discarding probes for which the associated gene was not represented in the RefSeq database 10 , and unifying all probes measuring the same gene (ignoring the potential differences among splice variants) yielded probes representing 7500 human genes. The experiments measuring replicates of the same biological condition were averaged to reduce the measurement noise, resulting in 47 data points per probe. We considered that a probe was expressed in a certain condition if its average reading was above a certain cutoff value. The results were not sensitive to the exact cut-off value, and we chose 200 standard Affymetrix averagedifference units, considered to be a conservative cut-off value for determining gene presence 9 . This is also the trimmed average expression level in each tissue in accordance with the standard Affymetrix normalization procedure 11 ; 12 . Thus, our HK genes are expressed in all tissues at an above-average level.

A histogram (Fig. 1) of the number of genes expressed in exactly M of the 47 tissues shows a clear tendency for frequency to decrease as M increases. However, a substantial number of genes (575), belong to the class of genes that are expressed in all tissues. Because their number is far greater than expected based on the general trend described above, we assumed this class to be rich in HK genes, and considered it to be the set of HK genes.

It is noteworthy that the genes in our HK list tend to have an average expression significantly higher than other genes; the geometric mean expression of our HK genes is 1200 in Affymetrix average difference units, whereas that of other genes is 150. The difference cannot be accounted for by the cutoff used to define the HK genes, and is not a result of a bias due to inclusion of genes expressed in a few tissues only (data not shown).

Two additional tests were conducted to validate this set. First, a study of the GO annotation 13  of these genes revealed the set is rich in metabolic proteins (24%) and RNA-interacting proteins (19%, mostly ribosomal proteins). Second, we compiled a list of 18 well-established HK genes commonly used for quantitative PCR calibration 14 ; 15 , and checked our list against it. We found 13 of the 18 genes in our list, and the other five were not represented on the microarray (see Table in Supplementary Information at http://www.compugen.co.il/supp\_info/Housekeeping\_genes.html).

## II Length analysis of HK genes

Table 1 compares the lengths of various parts of the HK genes and the background genes. The alignment data was taken from the UCSC genome browser (http://genome.ucsc.edu) 16 . We excluded 322 genes that do not have a unique alignment, as well as 1242 genes that were not expressed in any tissue (to avoid potential problems because of defective probes). This left 532 HK genes and 5404 non-HK genes. The histograms in Fig. 2-4 compare HK genes with the other genes by total intron length, 5’ UTR length and coding sequence length. Remarkably, there was a statistically significant difference between HK and non-HK genes in all aspects of gene structure. Average intron length is shorter for the HK genes than for the background genes (2573 bp versus 5025 bp, respectively); total gene length is shorter (21,050 bp versus 53,418 bp); average exon length is shorter (212 bp versus 240 bp); average lengths of both 3’ and 5’ untranslated regions (UTRs) are shorter (5’: 135 bp versus 173 bp; 3’: 599 bp versus 846 bp); and, most notably, the translated proteins are shorter as well (403 amino acids versus 590 amino acids). Accordingly, the number of introns bp per unit of coding sequence length is lower for the HK genes (20 versus 32). We studied the structure of each gene as a function of the number of tissues it is expressed in and verified that the results are not due to bias of the non-HK genes by tissue-specific genes (data not shown).

The pronounced statistical characteristics of the HK gene set further supports their assignment as a unique set. Our findings confirm and extend previous research, showing that the introns of highly expressed genes are shorter 5 . As mentioned above, the HK genes expression levels are high, and the fact that they have to be expressed in all cells at all times makes them even more costly to transcribe. Previously 8 , the high abundance of a certain gene in EST libraries was an indication the gene was highly expressed in the human body. It was pointed out 8 , however, that this method is prone to bias due to the inclusion of normalized and tumor libraries and overrepresentation of certain tissues. Our approach overcomes this difficulty and confirms the previous result. Moreover, we find here that UTRs and even the encoded proteins are shorter for the HK genes. The magnitude of the difference is greater for the introns than for the exons and proteins (Table 1), which makes sense because the coding sequences and the UTRs are less susceptible to change.

It should be mentioned that intronless genes were included in our analysis after verifying that their inclusion or exclusion had no effect on the results. It also must be noted that the UTRs are not always fully sequenced, and thus their actual lengths might be longer. This bias was found to have no effect on the length of the coding sequences, and in any case the effect would be the same for both HK and non-HK genes.

It has been noted that codon usage bias in nonmammalian organisms is correlated with the expression level and with the gene length 17 ; 18 ; 19 . These results led to the conjecture of selective pressure on highly expressed genes resulting in shorter proteins 19 . However, no evidence for this selection was found 18 , possibly because of a lack of high quality databases for these organisms. Recent works have suggested that there is no selection for codon usage bias in humans 20 , and thus our results demonstrate that the expression-length correlation is not related to the expression-codon bias correlation.

It could be argued that selection towards shorter genes should have eliminated the introns in highly expressed genes. However, it is known that introns do have important roles, such as splicing regulation. Therefore, there is a balance between the advantageous contribution of the introns and the selective pressure for shortening.

Finally, when we compared our results with two (largely overlapping) published sets of HK genes, we found that roughly half of the genes in the intersection of those sets were present in our set. We used the genomic structure to test the remaining genes, and found a statistically significant difference between them and our HK gene set. The differences between our results and those of earlier studies 3 ; 4  could be due to the fact that the database we used was based on more advanced chip technology and included many more different tissues, giving it more discriminative power to identify HK genes.

In conclusion, we have identified a set of HK genes. The set is publicly available at http://www.compugen.co.il/supp\_info/Housekeeping\_genes.html and can be used for calibration of microarrays, toxicity evaluation and quantitative PCR experiments. Furthermore, we show that HK genes have shorter introns, UTRs and coding sequences, attesting to the strong selection for compactness in these genes.

###### Acknowledgements.

We thank Andrew Su for helpful discussion and for providing us with the RefSeq mapping. Gady Cojocaru and Rotem Sorek are acknowledged for comments on the manuscript and insightful discussion.


text:
# Auxiliary fields and hadron dynamics

## I Introduction

A system of relativistic particles can be quantized as a constrained system by using the auxiliary field formalism (also known as einbein field formalism) to get rid of the square root term in the Lagrangian polyakov ; brin77 . Applied to the QCD Lagrangian, this technique yields Hamiltonians for mesons or baryons $`H(\mu _i)`$ depending on auxiliary fields $`\mu _i`$ simonov ; kala97 ; naro02 . The way of removing these fields for physical applications is not well defined and can lead to ambiguities. In practice, these fields, representing for instance the energy density for the particles, are finally treated as $`c`$-number and determined from the minimum energy conditions

$$\frac{E(\mu _i)}{\mu _i}=0,$$
(1)
where $`E(\mu _i)`$ is an eigenvalue of $`H(\mu _i)`$.

When the interaction $`V`$ between particles does not depend on $`\mu _i`$, the values of $`\mu _i`$ which minimize the hamiltonian operator $`H(\mu _i)`$ (and not $`E(\mu _i)`$) are $`\mu _i=\sqrt{\stackrel{}{p}_i^2+m_i^2}`$. The eigenvalue equation reduces then to a spinless Salpeter equation with the potential $`V`$. So, in this case, the auxiliary field formalism can be considered as an approximation of the spinless Salpeter equation luch96 .

Although handling a relativistic expression for the kinetic energy term is not a serious problem with nowadays numerical algorithms, a nonrelativistic expression is always simpler, specially for many body calculations. In a sense the auxiliary field formalism simulates a relativistic expression with the simplicity of a nonrelativistic one. The price to pay is a further minimization on eigenenergies. The values of $`\mu _i`$ are generally different for various radial excited states and this leads to non-orthogonality among them. This is another price to pay. An aim of this paper is to study whether these drawbacks are very inconvenient or not.

In Sec. II, we consider the case of two identical particles interacting via a potential independent of auxiliary fields. We solve the eigenvalue equations for hamiltonians written in the auxiliary field formalism and for corresponding hamiltonians in spinless Salpeter equations. Then, we compare the eigenvalues and eigenvectors for these two approaches and discuss the problem of non-orthogonality. The study is first made for a toy model where most of the results are analytical and then in a more realistic one which relies on a precise numerical treatment.

It is now believed that the color field between a quark and an antiquark in a meson can be approximated by a rigid string, carrying both energy and angular momentum. Some years ago, a phenomenological model taking into account this particular string dynamics has been developed: the relativistic flux tube model laco89 ; olss95 . The classical equations of motion take the form of a coupled system of nonlinear equations in which a particular quantity is introduced, the transverse velocity of the ends of the string, where the particles are located. More recently, a lagrangian for a rotating string was derived from the QCD lagrangian dubi94 . Adopting a simple straight line ansatz for the string, another coupled system of nonlinear equations was obtained, containing auxiliary functions for the particle energy density and the string energy density morg99 . In this case, the potential does depend on the auxiliary fields. In Sec. III, we show that these two approaches are exactly the same, provided all auxiliary functions are correctly eliminated from the rotating string hamiltonian.

## II Auxiliary fields and Spinless Salpeter equation

In the formalism of auxiliary fields, the hamiltonian $`H^A`$ for two particles with the same mass $`m`$, interacting via a potential $`V`$ which does not depend on $`\mu `$, is written kala97  ($`\mathrm{}=c=1`$)

$$H^A=\frac{\stackrel{}{p}^2+m^2}{\mu }+\mu +V.$$
(2)
The philosophy of this formalism is to keep the simplicity of the nonrelativistic form of the kinetic energy operator, to get the eigenenergies in terms of auxiliary fields, and to get rid of them through a minimization procedure.

As mentioned, the eigenvalues and eigenvectors of this hamiltonian depends on $`\mu `$

$$H^A|A;nl;\mu =E_{nl}^A(\mu )|A;nl;\mu ,$$
(3)
where $`l`$ is the orbital angular momentum and $`n`$ (0, 1, …) the vibrational quantum number. The physical solutions are obtained by choosing for each eigenstate the value of $`\mu `$ which minimizes its energy. We will note $`\mu _{nl}`$ the value of $`\mu `$ minimizing the eigenvalue $`E_{nl}^A(\mu )`$. The following notations will be used

$$E_{nl}^A=E_{nl}^A(\mu _{nl})E_{nl}^A(\mu )\mu \text{and}|A;nl=|A;nl;\mu _{nl}.$$
(4)
Applying the Hellmann-Feynman theorem feyn39  to $`H^A`$, we find

$$\frac{E_{nl}^A(\mu )}{\mu }=\frac{A;nl;\mu |\stackrel{}{p}^2+m^2|A;nl;\mu }{\mu ^2}+1.$$
(5)
It follows immediately that

$$\mu _{nl}^2=A;nl|\stackrel{}{p}^2+m^2|A;nl,$$
(6)
and that

$$E_{nl}^A=2\mu _{nl}+A;nl|V|A;nl.$$
(7)
Applied to $`H^A`$, the virial theorem yields $`2\stackrel{}{p}^2/\mu =\stackrel{}{r}\stackrel{}{}V`$. If the potential V is a homogeneous function of $`r`$ of degree $`\alpha `$, the formula (7) can be simplified. We have then

$$E_{nl}^A=\frac{2(1+\alpha )}{\alpha }\mu _{nl}\frac{2m^2}{\alpha \mu _{nl}}\text{if}\stackrel{}{r}\stackrel{}{}V=\alpha V.$$
(8)
The scalar product $`A;n^{}l^{}|A;nl`$ is proportional to $`\delta _{l^{}l}`$ because of the orthogonality of the spherical harmonics. But, for a fixed value of $`l`$, two eigenstates with different values of $`n`$ are not orthogonal since they are characterized by different values of $`\mu `$. So, it is interesting to compute the quantity which measures the overlap of two states

$$P_{n^{}nl}=|A;n^{}l|A;nl|^2.$$
(9)
We have obviously $`P_{nn^{}l}=P_{n^{}nl}`$ and $`P_{nnl}=1`$, and we can expect that $`P_{n^{}nl}\delta _{n^{}n}`$.

Another way of dealing with auxiliary field is to minimize not the eigenvalues but the hamiltonian operator itself and impose the condition $`H^A/\mu =0`$, $`\mu `$ being considered as an operator by its own. In the case of hamiltonian (2), it is easy to determine that the value of the auxiliary field $`\mu `$ which realizes this condition is $`\mu =\sqrt{\stackrel{}{p}^2+m^2}`$. The minimal hamiltonian (2) then becomes

$$H^S=2\sqrt{\stackrel{}{p}^2+m^2}+V,$$
(10)
which is simply a spinless Salpeter hamiltonian with the same potential $`V`$. We note its eigenvalues and eigenstates by

$$H^S|S;nl=E_{nl}^S|S;nl.$$
(11)
It has been shown in Ref. luch96 , that we always have $`E_{nl}^AE_{nl}^S`$.

If the two ways of treating the auxiliary field were expected to be very similar, one sees that the first method simulates the effect of a relativistic kinetic energy term with a nonrelativistic one; the potential is kept unchanged. The differences $`E^AE^S`$ give an indication on the quality of this procedure.

Since the minimal energy eigenvalues of $`H^A`$ are upper bound of the eigenvalues of $`H^S`$ and since these two hamiltonians have in some sense the same physical content, their corresponding eigenstates must be very similar. By computing the overlap

$$T_{nl}=|A;nl|S;nl|^2,$$
(12)
we can expect that $`T_{nl}1`$.

If we define

$$M_{nl}=S;nl|\sqrt{\stackrel{}{p}^2+m^2}|S;nl,$$
(13)
we can also expect that $`M_{nl}\mu _{nl}`$, because of the relation (6). In both models, $`m`$ is sometimes considered as the current mass, and $`M_{nl}`$ or $`\mu _{nl}`$ are then interpreted as the constituent mass, which, in this formalism, is state dependent.

For two massless particles, the relativistic virial theorem luch90  yields $`2\sqrt{\stackrel{}{p}^2}=\stackrel{}{r}\stackrel{}{}V`$. If the potential V is a homogeneous function of $`r`$ of degree $`\alpha `$, new relations can be written. We have then

$$E_{nl}^S=(1+\alpha )S;nl|V|S;nl=\frac{2(1+\alpha )}{\alpha }M_{nl}\text{if}\stackrel{}{r}\stackrel{}{}V=\alpha V\text{and}m=0.$$
(14)
When the mass $`m`$ increases, the differences existing between the two approaches $`H^A`$ and $`H^S`$ tend to vanish and $`\mu _{nl}m`$ for all values of quantum numbers. In order to quantify these differences, we will study two models containing a confining potential, which is the relevant kind of interaction for the hadronic physics. In particular, we will focus our attention on the case $`m=0`$, for which we expect the largest discrepancies between the two approaches.

### II.1 Toy model

We first consider a “toy model” for which most of the calculations are analytical. Let us choose $`V=kr^2`$ and $`m=0`$. In this case, $`H^A`$ is a harmonic oscillator hamiltonian. Its solutions are well-known and we have

$$E_{nl}^A(\mu )=2\sqrt{\frac{k}{\mu }}(2n+l+3/2)+\mu .$$
(15)
The regularized radial part of the wave function is noted

$$u_{nl}^A(r)=(\mu k)^{3/8}rO_{nl}\left((\mu k)^{1/4}r\right),$$
(16)
where $`O_{nl}(x)`$ is the usual radial harmonic oscillator wave function and where dimensioned factors ensure that the eigenstates $`|A;nl;\mu `$ are normalized.

The values of the parameter $`\mu `$ minimizing the energies are given by

$$\mu _{nl}=k^{1/3}(2n+l+3/2)^{2/3},$$
(17)
and the calculation of the eigenenergies yields

$$E_{nl}^A=3k^{1/3}(2n+l+3/2)^{2/3}=3\mu _{nl}.$$
(18)
The relation $`E_{nl}^A=3\mu _{nl}`$ can also be obtained directly from formula (8).

The overlap of two eigenstates can be computed easily. Knowing the values of $`\mu _{nl}`$, we have

$$P_{n^{}nl}=F_{n^{}nl}^2\left(\left(\frac{2n+l+3/2}{2n^{}+l+3/2}\right)^{1/6}\right),$$
(19)
where the function $`F_{n^{}nl}(x)`$ is given, with its properties, in Ref. sema95a . The quantity $`P_{n^{}nl}`$ is independent of $`k`$ and decreases with increasing values of $`l`$ and $`|nn^{}|`$. For instance, we have $`P_{\mathrm{1\hspace{0.17em}0\hspace{0.17em}0}}0.0287`$, which is the worse case. Considering the states as orthogonal is thus always a good approximation.

In order to solve the corresponding hamiltonian $`H^S`$ for our toy model, let us consider the following nonrelativistic hamiltonian

$$H^N=\frac{\stackrel{}{p}^2}{2\nu }+ar.$$
(20)
Solutions for $`l=0`$ are known luch91 . The eigenenergies are given by

$$E_{n0}^N=\left(\frac{a^2}{2\nu }\right)^{1/3}x_n,$$
(21)
where $`x_n`$ is the $`(n+1)`$th zero of the Airy function $`\text{Ai}(x)`$ ($`x_02.338`$, $`x_14.088`$, …). The regularized and normalized radial part of the eigenvectors is written as brau98

$$u_{n0}^N(r)=(2\nu a)^{1/6}A_n\left((2\nu a)^{1/3}r\right)\text{where}A_n(x)=\frac{\text{Ai}(x+x_n)}{\sqrt{_{x_n}^{\mathrm{}}\text{Ai}^2(q)𝑑q}}.$$
(22)
Let us note that we have $`_{x_n}^{\mathrm{}}\text{Ai}^2(q)𝑑q=\text{Ai}_{}^{}{}_{}{}^{2}(x_n)`$.

The hamiltonian $`H^N`$ is converted into the hamiltonian $`H^S`$ with $`V=kr^2`$ and $`m=0`$ by the duality transformation $`|\stackrel{}{p}|{\displaystyle \frac{1}{2}}ar`$, provided the parameters $`k`$, $`a`$, and $`\nu `$ are related by $`k={\displaystyle \frac{a^2}{8\nu }}`$ luch91 . Due to this relation, we have immediately the property

$$E_{n0}^S=(4k)^{1/3}x_n.$$
(23)
In order to compare energies $`E_{n0}^A`$ and $`E_{n0}^S`$ given respectively by relations (18) and (23), we can use the following approximation for the numbers $`x_n`$ abra70

$$x_n\left[\frac{3\pi }{2}(n+3/4)\right]^{2/3}.$$
(24)
With this formula, whose error decreases when $`n`$ increases, relation (23) can be written

$$E_{n0}^S\left(\frac{3\pi }{2}\right)^{2/3}k^{1/3}(2n+3/2)^{2/3}.$$
(25)
So we have immediately the ratio (independent of $`k`$),

$$\frac{E_{n0}^S}{E_{n0}^A}\frac{1}{3}\left(\frac{3\pi }{2}\right)^{2/3}0.937.$$
(26)
Numerically, we found that the ratio $`E_{n0}^S/E_{n0}^A`$ is around 0.944 for $`n=0`$ and tends rapidly towards the asymptotical value (26) as $`n`$ increases.

From relation (22) and the duality transformation, the regularized function of S-wave eigenvectors of $`H^S`$ can be obtained in the momentum space

$$u_{n0}^S(p)=\left(\frac{2}{k}\right)^{1/6}A_n\left(\left(\frac{2}{k}\right)^{1/3}p\right).$$
(27)
Consequently, the quantity $`M_{n0}=S;nl|\sqrt{\stackrel{}{p}^2}|S;nl`$ can be calculated by the following integral ($`\sqrt{\stackrel{}{p}^2}=|\stackrel{}{p}|`$)

$$M_{n0}=k^{1/3}\frac{1}{2^{1/3}}_0^{\mathrm{}}𝑑qqA_n^2(q)=(4k)^{1/3}\frac{x_n}{3}=\frac{E_{n0}^S}{3},$$
(28)
and we find that $`M_{n0}/\mu _{n0}=E_{n0}^S/E_{n0}^A0.937`$. The relation $`E_{n0}^S=3M_{n0}`$ can also be obtained directly from formula (14).

Since the Fourier transform of a harmonic oscillator is also a harmonic oscillator (with a phase factor), the overlap $`T_{n0}`$ can be computed in the momentum space. We obtain

$$T_{n0}=\frac{2^{1/3}}{\sqrt{2n+3/2}}\left[_0^{\mathrm{}}𝑑qqO_{n0}\left(\frac{q}{(2n+3/2)^{1/6}}\right)A_n\left(2^{1/3}q\right)\right]^2,$$
(29)
which is independent of $`k`$. We have checked numerically that $`T_{\mathrm{0\hspace{0.17em}0}}0.997`$ and that $`T_{n0}`$ decreases monotonically when $`n`$ increases. For instance, $`T_{\mathrm{9\hspace{0.17em}0}}0.523`$. One sees that, when $`n`$ increases, the eigenvalues for the two formalisms are almost the same (to within 6%), whereas the corresponding wave functions can differ appreciably.

### II.2 Realistic model

The previous toy model was interesting because most of the quantities can be calculated analytically. However, it lacks some physics. First the real confining potential is more alike a linear one, and second there exist one gluon exchange contributions which add a coulombic term, important at short distance. In order to stick more to the true physical situation, we switch, in this section, to a more realistic model. Of course, this study needs a complete numerical treatment.

In the framework of a semirelativistic potential model, a more realistic interaction to simulate the dynamics between a quark and an antiquark inside a meson is certainly what is called the funnel potential fulc94

$$V(r)=\frac{\kappa }{r}+ar\text{with}\kappa =0.5\text{and}a=0.2\text{GeV}^2.$$
(30)
The values chosen here for the parameters $`\kappa `$ and $`a`$ can be considered as typical. The eigenvalue equations for hamiltonians $`H^A`$ and $`H^S`$ with this potential have been numerically solved by the Lagrange-mesh method sema01 . This technique is very accurate and can be implemented easily for both nonrelativistic or semirelativistic kinematics.

In Fig. 1, the quantities $`\mathrm{\Delta }E=(E^AE^S)/E^S`$ are presented, as a function of the mass $`m`$, for the three lowest states (1S, 2S, 1P) obtained with the potential (30). As expected, $`\mathrm{\Delta }E`$ decreases when the quark mass increases. The corresponding curves for the quantities $`E^AE^S`$ present the same profiles. Thus we only present more results for the case of a vanishing quark mass which maximizes the error.

In Table 1, some energies $`E_{nl}^S`$ and $`E_{nl}^A`$ are given for $`n3`$ and $`l2`$, for the potential (30) and $`m=0`$. The difference $`E_{nl}^AE_{nl}^S`$ can be quite large, around 100 MeV for the ground state. It increases for increasing values of $`n`$ and decreases for increasing values of $`l`$. The quantity $`M_{nl}`$ is also given in this Table; it increases with $`n`$ and $`l`$. The parameter $`\mu _{nl}`$ can differ from $`M_{nl}`$ by several tens of MeV, but the curves $`E_{nl}^A(\mu )`$ always present a flat minimum around $`\mu _{nl}`$. So it is not necessary to obtain a precise value of $`\mu _{nl}`$ to obtain a good value of $`E_{nl}^A`$.

The overlap $`P_{nn^{}l}`$ is given in Table 2 for $`n,n^{}3`$ and for $`l1`$, for the potential (30) and $`m=0`$. Except in the case $`|nn^{}|=1`$, the overlap is very small. So the eigenstates $`|A;nl`$ can be considered as quasi-orthogonal. This conclusion is the same than in the toy model, and we can hope that it is a universal conclusion, valid whatever the potential $`V`$.

In table 3, the overlap $`T_{nl}`$ is presented for $`n3`$ and $`l2`$, for the potential (30) and $`m=0`$. It increases for increasing values of $`l`$ and decreases for increasing values of $`n`$, and present similar features to the toy model. For small values of the vibrational quantum number, the eigenvectors $`|A;nl`$ and $`|S;nl`$ are rather similar.

## III Rotating QCD string

Up to now, the confining term was put by hand. There exist more serious explanations of this term, based on QCD arguments. In particular, one can consider a string with an energy density between the quark and the antiquark. This string is itself a dynamical object. Here we examine two such approaches.

We first quickly present the relativistic flux tube model, which is essentially phenomenological, and then the rotating string model within the formalism of auxiliary fields, which has a more firm foundation on QCD theory. Lastly, we show that these two models are completely equivalent.

### III.1 Phenomenological model

In the simplest version of the relativistic flux tube model laco89 , a quark and an antiquark, with the same mass $`m`$, move being attached with a rigid flux tube, assumed to be linear with a uniform constant energy density $`a`$. The system rotates in a plane with a constant angular velocity around the center of mass, which is assumed to be stationary. Denoting $`r`$ the distance between the two particles, $`p_r`$ the radial momentum $`\left(p_r^2={\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{r^2}}r\right)`$, $`L`$ the total angular momentum, $`v_{}`$ the transverse velocity of the quark or the antiquark (relative to the string direction), $`\gamma _{}=(1v_{}^2)^{1/2}`$, and $`U(r)`$ a potential used to describe dynamical effects coming from mechanisms other than the flux tube, the classical equations of motion of the system are given by laco89

$`{\displaystyle \frac{L}{r}}`$ $`=`$ $`v_{}\gamma _{}\sqrt{p_r^2+m^2}+arf(v_{}),`$ (31)
$`H`$ $`=`$ $`2\gamma _{}\sqrt{p_r^2+m^2}+ar{\displaystyle \frac{\mathrm{arcsin}v_{}}{v_{}}}+U(r),`$ (32)
where

$$f(x)=\frac{1}{4x^2}\left(\mathrm{arcsin}xx\sqrt{1x^2}\right),$$
(33)
is a very important function in the formalism. These equations have been generalized in the case of asymmetrical systems olss95 . Numerical solutions of the quantized versions of these equations have been obtained laco89 ; sema94 ; olss95 ; sema95b  (for a practical calculation $`L`$ must be replaced by $`\sqrt{l(l+1)}`$). Equations (31)-(32) are coupled nonlinear equations. In practice, the value $`v_{}(m,a,L;p_r,r)`$ is extracted from Eq. (31) and injected into Eq. (32), giving a $`L`$ dependent hamiltonian $`H(m,a,L;p_r,r)`$, which is diagonalized afterwards.

When the relativistic flux tube hamiltonian is supplemented by appropriate potentials (Coulomb-like, instanton induced effects) and when it is assumed that each extremity of the flux tube can give a constant energy contribution, rather good meson spectra can be obtained sema94 ; sema95b .

### III.2 Auxiliary fields for a rotating string

Starting from the QCD lagrangian, a lagrangian for a meson can be derived taking into account the dynamical degrees of freedom of the string morg99 . For a system containing a quark and an antiquark with the same mass, with the hypothesis of a straight line configuration for the minimal string, and introducing an auxiliary field $`\mu `$ for the energy density of the quarks and another auxiliary field $`\nu `$ for the energy density of the string, the following classical hamiltonian can be obtained morg99

$`H={\displaystyle \frac{p_r^2+m^2}{\mu (\tau )}}+\mu (\tau )`$ $`+`$ $`{\displaystyle \frac{L^2/r^2}{\mu (\tau )+2{\displaystyle _0^1}\left(\beta {\displaystyle \frac{1}{2}}\right)^2\nu (\beta ,\tau )𝑑\beta }}`$ (34)
$`+`$ $`{\displaystyle \frac{a^2r^2}{2}}{\displaystyle _0^1}{\displaystyle \frac{d\beta }{\nu (\beta ,\tau )}}+{\displaystyle _0^1}{\displaystyle \frac{\nu (\beta ,\tau )}{2}}𝑑\beta +U(r),`$
where $`U(r)`$ is a potential describing dynamical effects coming from mechanisms other than the string. In this formula, $`\tau `$ is the common proper time of the two particles and $`\beta `$ is a coordinate along the string. Within this formalism, the two auxiliary functions $`\mu `$ and $`\nu `$ are to be varied and to be found from the minimum of $`H`$.

Let us first remark that, if we neglect the integral in the denominator of the term depending on $`L^2`$ or if we put $`L=0`$, then $`H`$ is minimum for $`\nu =ar`$, independent of $`\beta `$. In this case, because of the relation $`p_r^2+{\displaystyle \frac{L^2}{r^2}}=\stackrel{}{p}^2`$, $`H`$ reduces simply to $`H^A`$ with $`V(r)=ar+U(r)`$. The linear confining term appears naturally in the formalism. A further minimization on $`\mu `$, as demonstrated in second section, gives rise to $`H^S`$ with the same potential $`V(r)`$.

This result is interesting and simple, but indeed one can make an exact minimization on the field $`\nu `$. A straightforwards calculation shows that the minimization condition $`H/\nu =0`$ is fulfilled if the field $`\nu `$ is set to $`\nu _0`$ with

$$\nu _0(\beta )=\frac{ar}{\sqrt{14y^2\left(\beta \frac{1}{2}\right)^2}},$$
(35)
where $`y`$ is to be found from the transcendental equation

$$\frac{L}{ar^2}=f(y)+\frac{\mu y}{ar}.$$
(36)
The function $`f`$ defined previously by relation (33) appears again curiously (for a practical calculation $`L`$ must also be replaced by $`\sqrt{l(l+1)}`$). Using expression (35), one obtains from relation (34)

$$H=\frac{p_r^2+m^2}{\mu }+\mu (1+y^2)+ar\frac{\mathrm{arcsin}y}{y}+U(r).$$
(37)
WKB solutions for the system (36) and (37) are obtained in Ref. morg99  for $`m=0`$ and $`U(r)=0`$. Regge trajectories are computed in agreement with experimental data.

### III.3 Equivalence between the two string models

In principle the way of solving the system (36)-(37) is the following. One extracts $`y(a,L;r;\mu )`$ from Eq. (36) and inject it into Eq. (37) leading to $`H(m,a,L;p_r,r;\mu )`$. This hamiltonian is then diagonalized and, for each state, the optimal value of $`\mu `$ must be determined to get the physical eigenenergies. In view of the form for various expressions, this procedure seems hopeless. Fortunately this is not so and one can go one step further and perform the minimization on $`\mu `$.

The partial derivative of Eq. (36) with respect to $`\mu `$ gives the following relation

$$\frac{y^2}{ar}+\frac{y}{\mu }\left[\frac{\mu y}{ar}+\frac{1}{2y\sqrt{1y^2}}\frac{\mathrm{arcsin}y}{2y^2}\right]=0,$$
(38)
in which the expression $`y/\mu `$ appears. Now imposing $`H/\mu =0`$, we obtain from Eq. (37)

$$\frac{p_r^2+m^2}{\mu ^2}+1+y^2+\frac{y}{\mu }\left[2\mu y+\frac{ar}{y\sqrt{1y^2}}ar\frac{\mathrm{arcsin}y}{y^2}\right]=0.$$
(39)
Combining the two last formulas, we obtain

$$\mu =\sqrt{\frac{p_r^2+m^2}{1y^2}}.$$
(40)
Using this relation, it can be verified that Eqs. (38) and (39) are well defined for $`y[0,1]`$.

Equation (40) giving $`\mu (y)`$ and the presence of the function $`f`$ (33) in Eq. (36) are the clues for the identification of both methods. Replacing $`y`$ by $`v_{}`$ and $`\mu `$ by $`\gamma _{}\sqrt{p_r^2+m^2}`$ transforms Eq. (36) into Eq. (31) and Eq. (37) into Eq. (32), proving thus the complete equivalence between both approaches. In passing we have now a clear idea of the physical content for the auxiliary field $`\mu `$ and of the mysterious variable $`y`$ appearing in the theory.

## IV Conclusion

Effective hamiltonians for hadron dynamics can be derived from the original QCD lagrangian. These hamiltonians can depend on auxiliary fields, representing for instance the quark energy density kala97 ; naro02 . When the interaction does not depend on these fields (exactly or by approximation), it is possible to reduce the effective hamiltonian to a semirelativistic one (spinless Salpeter).

In the first part of this paper, we have compared the eigenvalues and eigenvectors for hamiltonians written in the auxiliary field formalism $`H^A`$ and for the corresponding spinless Salpeter hamiltonians $`H^S`$. We have only considered the case of two particles with the same mass and two different confining interaction: A “toy” quadratic potential and a more realistic funnel potential. We have shown that the eigenvalues of $`H^A`$ are close to the eigenvalues of $`H^S`$, but the relative differences can exceed 10%. The overlap of corresponding eigenvectors for the two approaches is generally close to unity for ground states, but it deteriorates when the vibrational quantum number increases. It appears that, if precise calculations are desired, it is better to solve directly a spinless Salpeter equations: Accurate numerical techniques exist for two-body sema01 , and many-body suzu98 ; silv01  problems. Hamiltonians with auxiliary fields can be useful if one searches to obtain analytical results, as it is the case in Ref. naro02 . The non-orthogonality resulting in this formalism is not a serious problem, the radial excited states being quasi-orthogonal.

The phenomenological semirelativistic flux tube hamiltonian has been developed in order to take into account the string dynamics in a meson laco89 ; olss95 . Another model derived from the QCD lagrangian, the rotating string, has been more recently developed morg99 . But the rotating string hamiltonian depends on auxiliary fields.

In the second part of this paper, we have shown that the classical equations of motion of the rotating string reduce exactly to the classical equations of motion of the semirelativistic flux tube, provided all auxiliary fields are correctly eliminated. The complete equivalence of both models has two virtues: to provide a natural interpretation of the auxiliary field in term of the transverse velocity and to reinforce the relevance of the relativistic flux tube model (which was considered up to now as a phenomenological one), since it is equivalent to a model based on firm QCD grounds.

The set of coupled equations for these two models are very complicated to solve numerically laco89 ; sema94 ; olss95 ; sema95b . The difficulties are essentially the same for both approaches, with the supplementary necessity to perform a minimization within the rotating string formalism in order to determine the value of the auxiliary field for the quark energy density. It is then more interesting to work directly with the semirelativistic flux tube model.

## V Acknowledgments

The authors are grateful to NATO services which give us the possibility to visit ourselves and greatly facilitates this work, through the NATO grant PST.CLG.978710. B. Silvestre-Brac and C. Semay (FNRS Research Associate position) would like to thank the agreement CNRS/CGRI-FNRS for financial support. They also thank A. M. Badalian, F. Brau, Yu. S. Kalashnikova, and Yu. A. Simonov for useful discussions.