text:
# The manifold spectra and morphologies of EROsBased on observations made at the European Southern Observatory, Paranal, Chile (Program 70.A-0140(B)).

## 1 Introduction

Extremely Red Objects (EROs) are important probes of galaxy evolution. On one hand, they allow to constrain the formation epoch, the evolutionary pattern and the clustering of the oldest early-type galaxies at $`z`$ $`\stackrel{>}{}`$1 (e.g. Daddi et al. 2000a,b; McCarthy et al. 2001; Cimatti et al. 2002a; Smith et al. 2002; Firth et al. 2002; Moustakas & Somerville 2002; Roche et al. 2002, 2003; Saracco et al. 2003). On the other hand, they allow to select dusty star-forming galaxies and obscured AGN in a way complementary to submillimeter surveys (Cimatti et al. 1998; Afonso et al. 2001; Mohan et al. 2002; Cimatti et al. 2002a; Franceschini et al. 2002; Brusa et al. 2002). Although recent VLT spectroscopy shed new light on moderately bright EROs with $`K_\mathrm{s}`$$`\stackrel{<}{}`$19 (Cimatti et al. 2002a), the characteristics of fainter EROs is still unclear mainly because of their weakness which makes spectroscopy very challenging even with 8-10m-class telescopes. However, as different morphologies are expected in case of spheroids or star-forming systems, HST imaging can be used to complement spectroscopy for the faint ERO population (e.g. Moriondo et al. 2000; Yan & Thompson 2003). The main open questions about EROs include: the relative fraction of different ERO types, the nature and star formation rates of dusty EROs, the number density of the oldest high-$`z`$ spheroids selected as EROs, the fraction of obscured AGN hosted by EROs. In this paper we present the main results of very deep VLT spectroscopy, HST imaging and VLA observations aimed at investigating the nature of a complete sample of 47 EROs with $`K_\mathrm{s}<20`$. The properties of the sample will be presented in a companion paper (Cassata et al. 2003, in preparation). $`H_0=70`$ km s-1 Mpc-1, $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ are adopted.

## 2 The sample and the observations

The complete sample (already described by Cimatti et al. 2002a) is made by the 47 EROs with $`RK_\mathrm{s}>5`$ and $`K_\mathrm{s}<20`$ present in the 32.2 arcmin2 area of the Chandra Deep Field South (CDFS; Giacconi et al. 2001) covered by the K20 survey (Cimatti et al. 2002b). One object was removed from the original sample because it was found to be a blend of two objects both with $`K_\mathrm{s}>20`$.

New VLT-UT4+FORS2 multi-object optical spectroscopy was made in November 2002 with the red-sensitive CCD array. The total integration time was 7.8 hours during 0.7′′-1.0′′ seeing conditions, with 1.0′′ wide slits, grism 300I (R$``$660) and dithering of the targets along the slits. The mask included 13 EROs of the K20/CDFS sample that were previously unidentified or not observed. The spectra were reduced and analyzed with the methods described Cimatti et al. (2002a, 2002b). VLT spectroscopy was complemented with the HST+ACS imaging in the $`BViz`$ bands taken in the GOODS/HST Treasury Program (Giavalisco et al. 2003).

Deep radio data were obtained with the VLA at 1.4 GHz in 1999 October with a bandwidth of 25 MHz in each of two circular polarizations, and at 4.8 GHz in 1999 Feb with a bandwidth of 100 MHz in each of two circular polarizations, for 12-hours total at both frequencies (Kellermann et al. 2003, in preparation). The rms sensitivities were 8.6 $`\mu `$Jy and 9.6 $`\mu `$Jy at 1.4 and 4.8 GHz, respectively, each with a full-width half-power resolution of $`3.5^{\prime \prime }`$. The reductions were made using AIPS and each field was self-calibrated. The position and flux density of each source with a peak flux density greater than 3.5-rms was determined from a Gaussian model fit of the source.

## 3 VLT spectroscopy and HST imaging results

The 13 ERO spectra provided 9 new redshift measurements: 7 star-forming galaxies at $`z1.21.6`$ with \[OII\]$`\lambda `$3727 emission, and 2 early-type objects at $`z1.5`$. When added to the previous identifications of Cimatti et al. (2002a), the spectroscopic redshift completeness for the EROs in the K20/CDFS area increases to 81% (21/26) ($`K_\mathrm{s}<19.2`$) and 62% (29/47) ($`K_\mathrm{s}<20.0`$). Photometric redshifts for the 18 spectroscopically unidentified EROs were derived using the ESO/GOODS VLT+FORS1 $`BVRIz`$ and the new VLT+ISAAC $`JHK_\mathrm{s}`$ public imaging and the $`Hyperz`$ software (Bolzonella et al. 2000). The agreement between spectroscopic and photometric redshifts is good and does not depend on the ERO class, with $`\mathrm{\Delta }z=z_{\mathrm{spec}}z_{\mathrm{phot}}`$=0.001 and $`\sigma (\mathrm{\Delta }z)`$=0.16 (29 EROs). Figure 1 shows the redshift distributions divided into morphological classes (see Section 4).

The morphological type of EROs was first estimated by visual inspection of the ACS images and isophotal maps in the four available bands ($`BViz`$). One ERO (not observed spectroscopically) is unresolved and has colors consistent with being a M star. Since in some cases the spatial resolution and/or the S/N ratio turned out to be limited, we just assigned each object to one of the following broad classes: E/S0 galaxies (class 1), spiral-like objects (class 2) and irregular/interacting systems (class 3). Ambiguous cases were flagged by introducing the intermediate classes 1.5 and 2.5. Thus, we used the GASPHOT software (Pignatelli and Fasano 1999) to carry out a quantitative analysis on the $`z`$-band ACS images (filter F850LP) by fitting the surface brightness profile of each ERO with the Sersic law ($`SBr^{1/n}`$), convolved with the average PSF extracted from the stars in the field. Being the bulk of the sample at $`z1`$, the filter F850LP is close to the rest-frame $`B`$-band, thus minimizing the effects of dust extinction and morphological K-correction for the available set of ACS filters. EROs belonging to the class 1.5 and having $`n`$2.5 or $``$1.5 were assigned to the class 1 or 2, respectively, while the cases with 1.5$`<n<`$2.5 were assigned to intermediate class 1.5. Instead, the Sersic index showed to be ineffective in disentangling classes 2 and 3, making advisable to leave unchanged most of the dubious cases with flag 2.5 (10 objects). Indeed, for six of these objects the fitting algorithm failed to converge, due to the extremely irregular and noisy image structure. These objects are not reported in the lower panel of Fig.2, where the final classification is plotted versus the Sersic index. We note that the larger number of class-2.5 vs. class-1.5 objects and the higher redshifts of “irregular” EROs (Fig.1) may be due to morphological K-correction effects (e.g. Marcum et al. 2001; Windhorst et al. 2002). By taking into account the previous uncertainties, the ranges of morphological fractions are 30–37% (class 1 – class 1+1.5), 24–46% (class 2 – class 2+2.5) and 17–39% (class 3 – class 2.5+3). The fractions are in reasonable agreement with the recent results by Moustakas et al. (2003) and Gilbank et al. (2003), and do not significantly change as a function of the limiting $`K_\mathrm{s}`$ magnitude. Fig.2 (top panel) and Fig.3 show that a rather good agreement is present between the morphological and spectroscopic classes.

## 4 Dusty star-forming EROs and obscured AGN

Only one object was detected at $`>5\sigma `$ level ($`S>50`$ $`\mu `$Jy) in the 1.4 GHz radio source catalog of the CDFS (Kellermann et al. 2003, in preparation). Our radio detection rate (1$`{}_{0.83}{}^{}{}_{}{}^{+2.30}`$/46 at $`1\sigma `$ poissonian confidence level) is compatible with the preliminary results of the Phoenix Deep Survey (Hopkins et al. 2003) where $``$6% of a sample of $``$425 EROs ($`K_\mathrm{s}<20`$, $`RK_\mathrm{s}>5`$) are detected with fluxes $`>60`$ $`\mu `$Jy ($`5\sigma `$) (J. Afonso, private communication), with the results of Smail et al. (2002) (3/50 EROs with $`K_\mathrm{s}<20`$, $`RK_\mathrm{s}>5.3`$ and $`S>50`$ $`\mu `$Jy at 1.4 GHz) and Roche et al. (2002) (1/31 EROs with $`K_\mathrm{s}<19.5`$, $`RK_\mathrm{s}>5`$ and $`S>50`$ $`\mu `$Jy at 1.4 GHz).

The ERO detected in the radio is an emission line object with \[OII\]$`\lambda `$3727 at $`z=1.610`$ and a spiral-like morphology. The flux densities at 1.4 GHz and 4.8 GHz are 106.8$`\pm `$8.6 $`\mu `$Jy and 155.1$`\pm `$9.6 $`\mu `$Jy respectively, implying a flat/inverted spectrum with $`\alpha _{\mathrm{radio}}=+0.3`$ ($`S(\nu )\nu ^\alpha `$). The object is also present in the X-ray catalog of Alexander et al. (2003) (RA=03:32:10.80, Dec= -27:46:27.6, J2000) with fluxes F(0.5-2 keV)=8.67$`\times 10^{17}`$ cgs and F(0.5-8 keV)=3.08$`\times 10^{16}`$ cgs, corresponding to luminosities of 1.4$`\times 10^{42}`$ erg s-1 and 5$`\times 10^{42}`$ erg s-1 respectively. Its X-ray luminosity, $`\alpha _{\mathrm{radio}}`$ and location in typical segregation diagrams (e.g. $`\alpha _{\mathrm{radio}}`$ vs. radio flux, optical flux vs. 0.5-8 keV flux; Bauer et al. 2002) are more suggestive of a Type 2 AGN than a starburst galaxy. The lack of broad MgII$`\lambda `$2800 emission implies obscuration of the active nucleus. From the X-ray hardness ratio we derive log(NH)$`<`$ 22.7 cm-2 (assuming a photon index $`\mathrm{\Gamma }=1.7`$).

Although EROs are often found as faint counterparts of X-ray selected sources (e.g. Yan et al. 2003), the fraction of $`K`$-selected EROs with X-ray emission is rather small (e.g. Alexander et al. 2002; 14$`{}_{7}{}^{}{}_{}{}^{+11}`$% at $`K<`$20.1; Roche et al. 2003; $``$8% at $`K_\mathrm{s}<20`$). Taking into account the only other X-ray source detected in our ERO sample (Brusa et al. 2002), the fraction of X-ray emitting EROs amounts to 2/46 (1.5-10% at 1$`\sigma `$ poissonian confidence level), consistently with the above mentioned results.

Although only one ERO has been detected significantly in the radio, a mean flux density of the undetected sources was estimated by comparing the flux density distribution of the non-detected ERO sources with that expected from noise and confusion. The radio intensity distribution at 1.4 GHz for the 28 non-detected spiral-like + irregular ERO sources has a rms of $`9.3\mu `$Jy and an average value of $`3.5\mu `$Jy compared with an rms of $`8.9\mu `$Jy and average value of $`0.8\mu `$Jy for the control sample of 270 random locations in the field. From this excess we infer that the average flux density of ERO sources at 1.4 GHz is $`3.5\pm 1.7\mu `$Jy. At 5 GHz we estimate an average flux density of $`<4\mu `$Jy. For the 17 E/S0 galaxies we obtained limit of $`<3\mu `$Jy at both 1.4 and 5 GHz.

Adopting a radio spectral index $`\alpha =0.6`$, the average rest-frame 1.4 GHz luminosity is $`2.7\times 10^{22}`$ W Hz-1 for the average redshift of the 28 undetected star-forming EROs ($`z=1.337`$). Following Condon (1992) and Haarsma et al. (2000) (with Salpeter IMF and 0.1-100 M mass range), we derive SFR$``$ 33 M yr-1. The estimated SFRs are consistent with those obtained by Brusa et al. (2002) using the X-ray emission (5–44 M yr-1). Following Kennicutt (1998), the expected average far-IR luminosity is $`12\times 10^{11}`$ L corresponding to 850$`\mu `$m observed fluxes below 1 mJy, thus explaining the low detection rate of $`K_\mathrm{s}<20`$ EROs with submm/mm continuum observations (e.g. Mohan et al. 2002). In comparison, the \[OII\]$`\lambda `$3727 luminosity of all the spectroscopically identified star-forming EROs in the CDFS provides an average SFR=2.8 M yr-1. If this difference with the radio is solely due to dust extinction, the required average reddening is $`E(BV)0.5`$. This amount of reddening is significantly larger than in Lyman-break galaxies (median $`E(BV)0.17`$, range $`E(BV)0.10.3`$; Giavalisco 2002 and references therein). The average SFR(radio) implies a cosmic SFR density of EROs at $`z1.3`$ of $`>`$0.008 M yr-1 Mpc-3 (adopting the comoving volume between $`z_{\mathrm{min}}=0.726`$ and $`z_{\mathrm{max}}=2.185`$), corresponding to about $`>`$12% of the total budget at that redshift (e.g. Somerville et al. 2001). This is strictly a lower limit because our sample is limited to $`K_\mathrm{s}<20`$ and due to the conservative volume adopted between $`z_{\mathrm{min}}`$ and $`z_{\mathrm{max}}`$.

## 5 The oldest and massive spheroids at $`z`$ $`\stackrel{>}{}`$1

The reddest spheroidals selected as EROs represent the oldest envelope of the E/S0 population at $`z`$ $`\stackrel{>}{}`$1 and are crucial to constrain the earliest epoch of massive system formation.

Taking into account the non-poissonian fluctuations due to clustering (Daddi et al. 2000a; Roche et al. 2002), the observed surface density of E/S0 galaxies (morphological class 1-1.5) is $`0.50_{0.16}^{+0.37}`$ arcmin-2 ($`1\sigma `$) ($`0.50_{0.26}^{+1.50}`$ at $`2\sigma `$). The Pure Luminosity Evolution (PLE) model of Daddi et al. (2000b) with a star formation $`e`$-folding time $`\tau =0.3`$ Gyr and a formation redshift $`z_\mathrm{f}2.5`$ predicts surface densities of $``$2 arcmin-2 and $``$0.9 arcmin-2 for the 2MASS local early-type (E+S0+ some contamination by Sa galaxies; Kochanek et al. 2001) and the local elliptical (E only; Marzke et al. 1994) luminosity functions respectively. Although the large uncertainties on the surface density prevent a stringent comparison with the PLE model, we note that the observed density of E/S0 galaxies at 0.8$`\stackrel{<}{}`$$`z`$$`\stackrel{<}{}`$1.6 is still formally consistent with the above PLE predictions within 1-2$`\sigma `$. The question of high-$`z`$ early-type galaxies is still controversial (e.g. Treu 2003 and references therein), and large area surveys are required to overcome the effects of the strong clustering. Figure 4 shows a significant spread in the formation redshifts of the identified E/S0 galaxies (2.2$`\stackrel{<}{}`$$`z_\mathrm{f}`$$`\stackrel{<}{}`$4.0 for $`\tau =0.3`$ Gyr). If $`\tau =0.1`$ Gyr, the range is 2$`\stackrel{<}{}`$$`z_\mathrm{f}`$$`\stackrel{<}{}`$3. From the lowest acceptable formation redshift ($`z_\mathrm{f}2.0`$) to $`z1`$ there is enough time to reach an age of $``$2.5 Gyr, compatible with that estimated from the average spectrum and from the 4000 Å continuum break amplitude of old EROs at $`z1`$ (Cimatti et al. 2002a), and within the range of the ages of $`z>1`$ spheroidal galaxies ($``$1.5- 3.5 Gyr) estimated with different spectral synthesis models (e.g. Nolan et al. 2001 and references therein). The rest-frame near-IR luminosities of old EROs in our sample are in the range of -24.5 and -27 (Pozzetti et al. 2003), and their stellar masses are $`10^{11}`$$`10^{12}`$ M (Fontana et al. 2003, in preparation).

The morphology, spectral properties, luminosities, inferred ages, formation redshifts and stellar masses imply the existence of a substantial population of old, passively evolving and fully assembled massive spheroids at $`z11.5`$ which requires that major episodes of massive galaxy formation occurred at $`z`$$`\stackrel{>}{}`$2. Daddi et al. (2003) uncovered a population of massive and clustered star-forming galaxies at $`z2`$ which may represent the long sought-for progenitors of massive spheroids.

## 6 Testing the photometric segregation of EROs

We used our sample to test the segregation of morphologically classified EROs in the $`RK_\mathrm{s}`$$`JK_\mathrm{s}`$ diagram proposed by Pozzetti & Mannucci (2000) (Fig. 5). The colors were derived from the new public ESO GOODS imaging taken with the VLT FORS1 ($`R`$) and ISAAC ($`JK_\mathrm{s}`$).

The relation used in Fig. 5 (dashed line) to define the boundary line between “old” and “dusty” EROs is $`JK_\mathrm{s}`$=0.34($`RK_\mathrm{s}`$)+0.09. This relation was derived specifically for the set of $`RJK_\mathrm{s}`$ filters of our photometry, and is valid for $`RK_\mathrm{s}>`$5.3 (Pozzetti & Mannucci 2000).

Two main results are evident from Fig.5: (1) the E/S0 galaxies with secure morphological classification (class=1) and most of star-forming EROs lie in their expected regions, (2) a few star-forming EROs (classes=2-3), however, are found in the region of E/S0 galaxies. If we consider only the EROs with reliable morphological classification, the contamination of the E/S0 region is 4/15 (27%), whereas 4/18 (22%) star-forming EROs do not fall in their expected locus. The contamination of the E/S0 region may be due to star-forming galaxies with a substantial underlying old stellar population which dominates the colors. Our results indicate that the color – color segregation is efficient and that the typical contamination amounts to $``$20-30% (see also Pierini et al. 2003).

###### Acknowledgements.

We thank R. Gilli for useful discussion, J. Afonso for information about the Phoenix survey and the referees for the constructive suggestions.


format: text

text:
# ITEP-LAT/2003-26 On the gauge dependence of the singlet and adjoint potentials

## 1 Introduction

The free energy of a static quark–antiquark pair – static quark potential – is of great importance for the understanding of the confinement–deconfinement transition, long distance properties of chromoplasma and heavy quark phenomenology at finite temperatures. One expects that the static quark potential depends on the color channel one chooses, e.g. singlet or adjoint, and that these potentials are gauge invariant . It is well known that the correlator of the Polyakov loop gives rise to the gauge invariant definition of the color averaged potential $`V_{av}(\stackrel{}{R})`$ . However, there is a problem with the gauge invariant definition of the singlet $`V_{sing}(\stackrel{}{R})`$ and adjoint $`V_{adj}(\stackrel{}{R})`$ potentials. Recently, it has been argued  that it is possible to overcome this problem and to arrive at a gauge invariant definition of the singlet and adjoint potentials using the so called dressed Polyakov lines. The dressing of the source may be viewed as a gauge transformation, and this approach is equivalent to the definition of the singlet and adjoint potentials in a certain gauge. It has been claimed  that one can use any unique gauge that is local in time, i.e. the gauge that does not change the spectrum of the transfer matrix (contrary to the case of the Lorentz/Landau gauge used in ). Therefore, to find $`V_{sing}(\stackrel{}{R})`$ one must calculate the average $`\text{Tr}\left(L(\stackrel{}{R})L^{}(\stackrel{}{0})\right)`$ in the chosen (time local) gauge, $`L(\stackrel{}{x})`$ being the Polyakov line defined in the next section. This approach has been used to calculate $`V_{sing}(\stackrel{}{R})`$ and $`V_{adj}(\stackrel{}{R})`$ in $`3D`$  and $`4d`$  $`SU(2)`$, and in $`4d`$ $`SU(3)`$  lattice gauge theories.

The main goal of this note is to study more closely the dependence of $`V_{sing}(\stackrel{}{R})`$ and $`V_{adj}(\stackrel{}{R})`$ on the choice of the gauge. To this purpose we calculated the gauge invariant correlator $`\text{Tr}L(\stackrel{}{R})\text{Tr}L^{}(\stackrel{}{0})`$ as well as the averages $`\text{Tr}\left(L(\stackrel{}{R})L^{}(\stackrel{}{0})\right)`$ in various local in time gauges : maximal tree axial gauge (AG), Coulomb gauge (CG) and generalized Coulomb gauge (GCG).

All our calculations were made in $`4d`$ $`SU(2)`$ lattice gauge theory in the finite volume with periodic boundary conditions. Main definitions and computational details are given in the next section. In section 3 we present our results and make conclusions.

## 2 Potentials, gauges and details of simulations

Let $`N_\tau `$ and $`N_\sigma `$ be the number of sites along the timelike and spacelike directions, respectively, and $`a`$ be the lattice spacing. We choose $`0x_iN_\sigma aa`$ ($`i=1,2,3`$) and $`0x_4N_\tau aa`$.

In the case of $`SU(2)`$ lattice gauge theory the color averaged potential $`V_{av}(\stackrel{}{R})`$ is given by

$$e^{V_{av}(\stackrel{}{R})/T}=\frac{1}{4}\text{Tr}L(\stackrel{}{R})\text{Tr}L^{}(\stackrel{}{0}).$$
(2.1)
where $`T=1/aN_\tau `$ is the temperature, $`L(\stackrel{}{x})`$ is the Polyakov line defined as $`L(\stackrel{}{x})=\frac{1}{2}\mathrm{Tr}_{x_0=0}^{N_\tau aa}U_{x_0,\stackrel{}{x},0}`$, $`U_{x\mu }`$ is the lattice field. Following  one can introduce the color singlet and adjoint potentials as follows

$`e^{V_{sing}(\stackrel{}{R})/T}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\text{Tr}\left(L(\stackrel{}{R})L^{}(\stackrel{}{0})\right);`$ (2.2)
$`{\displaystyle \frac{3}{4}}e^{V_{adj}(\stackrel{}{R})/T}`$ $`=`$ $`e^{V_{av}(\stackrel{}{R})/T}{\displaystyle \frac{1}{4}}e^{V_{sing}(\stackrel{}{R})/T}.`$ (2.3)
Eq.’s (2.2),(2.3) have meaning only when some gauge fixing condition is imposed. Without gauge fixing $`V_{sing}(\stackrel{}{R})=V_{adj}(\stackrel{}{R})=V_{av}(\stackrel{}{R})`$. In this work we employ the following local in time gauges.

* $`3D`$ maximal tree axial gauge.

In this case the gauge is fixed at every time slice independently as follows :

+ first we fix $`U_{x1}=1`$ for all $`x_1,x_2,x_3`$ apart from $`x_1=N_\sigma aa`$;
    + then we fix $`U_{x2}=1`$ for $`x_1=N_\sigma aa`$ and all $`x_2,x_3`$ apart from $`x_2=N_\sigma aa`$;
    + finally we fix $`U_{x3}=1`$ for $`x_1=x_2=N_\sigma aa`$ and all $`x_3`$ apart from $`x_3=N_\sigma aa`$.
* Generalized Coulomb gauge.

To fix this gauge one should find the maximum of the functional

$$F_U^\lambda (\mathrm{\Omega })=\underset{x}{}\frac{1}{2}\text{Tr}U_{x1}^\mathrm{\Omega }+\lambda \underset{x}{}\left[\frac{1}{2}\text{Tr}U_{x2}^\mathrm{\Omega }+\frac{1}{2}\text{Tr}U_{x3}^\mathrm{\Omega }\right],$$
(2.4)
with respect to gauge transformations $`\mathrm{\Omega }`$, where $`U_{x\mu }^\mathrm{\Omega }=\mathrm{\Omega }_xU_{x\mu }\mathrm{\Omega }_{x+\mu }^{}`$ and $`\lambda `$ is some parameter taking values between zero and unity. Evidently, at $`\lambda =1`$ one arrives at the standard Coulomb gauge.

Our main simulations were made on $`6\times 18^3`$ lattice at $`\beta =2.35`$. At this $`\beta `$ $`T/T_c=\sqrt{\sigma (\beta _c)/\sigma (2.35)}=0.76(1)`$, where $`\sigma `$ is string tension, and $`\beta _c=2.43`$ 111We took $`\sigma (2.35)=0.311(2)`$ from  and $`\sigma (2.43)`$ was estimated by interpolation of the data for string tension taken from the literature.. 5,000 configurations, separated by 100 combined sweeps consisting of one heat bath and two overrelaxation sweeps, were collected. To compare our results for the Coulomb gauge with those of ref.  we also made simulations on $`4\times 16^3`$ lattice at $`\beta =2.234`$ ($`T/T_c=0.80`$). To fix the Coulomb and the generalized Coulomb gauges we used the ordinary procedure alternating the relaxation and overrelaxation sweeps.

Statistical errors are calculated by the blocked jackknife method. To see the dependence of the potentials on the choice of the Gribov copy in the Coulomb gauge we generated $`5`$ random gauge copies for every equilibrium configuration at $`\beta =2.234`$. No sizable effect of the Gribov copies has been found.

## 3 Results in various gauges

In the Coulomb gauge our results are in agreement with earlier results  for the $`4d`$ case. The singlet potential $`V_{sing}^{CG}(\stackrel{}{R})`$ is below the color averaged $`V_{av}(\stackrel{}{R})`$ at small distances and approaches it from below with increasing $`|\stackrel{}{R}|`$, while the adjoint potential $`V_{adj}^{CG}(\stackrel{}{R})`$ is above $`V_{av}(\stackrel{}{R})`$ at small distances and approaches it from above at large distances (see Figure 1).

Of special interest for us is the comparison of the singlet and adjoint potentials in various time local gauges. Figure 1 shows the behavior of the singlet potential $`V_{sing}^{AG}(\stackrel{}{R})`$ in the maximal tree axial gauge (filled circles). The important observation is that this potential is far from being equal to the singlet potential $`V_{sing}^{CG}(\stackrel{}{R})`$ defined in the Coulomb gauge (triangles up Figure 1). In fact, this observation represents one of the main results of our work.

Numerically, the potential $`V_{sing}^{AG}(\stackrel{}{R})`$ is very close to the color averaged potential $`V_{sing}^{CG}(\stackrel{}{R})`$ for most values of the distance $`R`$. At the same time at distances along the lattice axes one can see clear deviations from $`V_{av}(\stackrel{}{R})`$ which are due to the lack of the rotational invariance of maximal tree axial gauge. The effects of this lack of the rotational invariance look even more impressive for the potentials measured along the axes in the $`1`$st, $`2`$nd and $`3`$rd directions separately (we do not show it in the Figure).

In this gauge the correlator eq.(2.2) can be explicitly expressed as a linear combination of the periodic Wilson loops (PWL), introduced in . The path connecting points $`\stackrel{}{0}`$ and $`\stackrel{}{R}`$ (space-like leg of PWL) depends on direction, e.g. the path is a straight line for $`\stackrel{}{R}=\{R,0,0\}`$, while it deviates from a straight line for $`\stackrel{}{R}=\{0,R,0\}`$ and $`\stackrel{}{R}=\{0,0,R\}`$. Thus the lack of the rotational invariance can be ascribed to the strong dependence of the PWL on this path. Let us also mention that the correlator eq.(2.2) in AG is exactly the same in the maximal tree gauge, which includes all four directions.

In the case of another gauge we used – GCG – the gauge fixing functional defined in eq.(2.4) depends on the parameter $`\lambda `$ chosen to be $`0<\lambda 1`$. Evidently, the choice $`\lambda =1`$ corresponds to the standard Coulomb gauge. We found that in this gauge both potentials, i.e the singlet potential $`V_{sing}^\lambda (\stackrel{}{R})`$ and adjoint potential $`V_{adj}^\lambda (\stackrel{}{R})`$, demonstrate a nontrivial dependence on $`\lambda `$.

In Figure 2 we show the $`R`$–dependence of the singlet potential $`V_{sing}^\lambda (\stackrel{}{R})`$ in GCG at some comparatively small values of $`\lambda `$ as well as $`R`$–dependence of $`V_{sing}^{CG}(\stackrel{}{R})`$ ($`\lambda =1`$) and $`V_{av}(\stackrel{}{R})`$. One can see that with decreasing $`\lambda `$ the singlet potential $`V_{sing}^\lambda (\stackrel{}{R})`$ moves towards the color averaged potential $`V_{av}(\stackrel{}{R})`$.

Figure 3 shows the dependence on $`\lambda `$ of the singlet potential $`V_{sing}^\lambda (\stackrel{}{R})`$ at some particular values of $`R`$ and reveals a nontrivial character of this dependence. Indeed, the data suggest that for $`\lambda \lambda _c`$ ($`\lambda _c0.6`$) there is no any sizable dependence on $`\lambda `$ and, therefore, $`V_{sing}^\lambda (\stackrel{}{R})=V_{sing}^{CG}(\stackrel{}{R})`$ within our errorbars. At $`\lambda <\lambda _c`$ this is not true anymore and the deviation of $`V_{sing}^\lambda (\stackrel{}{R})`$ from $`V_{sing}^{CG}(\stackrel{}{R})`$ increases with decreasing $`\lambda `$. Thus we find the class of the gauges, namely the GCG for $`\lambda \lambda _c`$, which give rise to the same singlet potential as the Coulomb gauge does. On the other hand, we find further evidence for the gauge dependence of the correlator eq.(2.2).

It is worthwhile to note that the ”transition” at $`\lambda =\lambda _c`$ has nothing to do with the restoration of the rotational invariance. Indeed, we found that at $`\lambda >\lambda _c`$ the singlet potential measured along the $`1`$st direction is different from the potentials measured along the $`2`$nd or $`3`$rd directions.

At the end we want to make a remark about the calculation of the singlet potential at $`T=0`$ using gauge noninvariant correlators of the Wilson lines. It was demonstrated in  that the singlet potential calculated from such correlators in the local in time Laplacian gauge was in full agreement with the gauge independent calculation using Wilson loops 222The correlator of the Wilson lines summed over spatial 2-plane was used to calculate the string tension in SU(3) theory in .. We have calculated the correlator of the Wilson lines in GCG at $`\lambda =0.2`$ and in CG on $`12^4`$ lattice at $`\beta =2.3`$. We found that the singlet potentials in both gauges agreed nicely with the potential determined from the Wilson loops. Thus at $`T=0`$ the gauge dependent correlators of the Wilson lines give rise to the same singlet potentials in GCG as in the Laplacian gauge or in CG. In AG the noise was too strong to make any definite conclusion.

We conclude that the statement made in  about uniqueness of the singlet and adjoint potentials defined by eq.’s (2.2),(2.3) in any unique and local in time gauge is possibly too strong. We believe that further study of the definition of the singlet and adjoint potentials at finite temperature will be useful. Our study can be extended to other, rotationally invariant gauges, different from the Laplacian gauge and the Coulomb gauge. An example of such gauge is the 3D Maximal Abelian gauge  with additional fixing of the abelian degrees of freedom or its analog using corresponding Laplacian operators.

This work has been supported by the grant INTAS-00-00111, RFBR grant 02-02-17308 and Heisenberg–Landau program.


text:
# Triggering star formation by galaxy-galaxy interactions

## References

Alonso, M. S, Tissera, P. B., Coldwell, G., Lambas, D. G., 2003, submitted.

Barton E. J., Geller M. J., Kenyon S. J., 2000, ApJ, 530, 660

Lambas, D. G., Tissera, P. B., Alonso, M. S. Coldwell, G. 2003, MNRAS in press


text:
# Evolutionary sequences of rotating protoneutron stars

## 1 Introduction

The first minute of life of non-rotating protoneutron stars (PNS) has been studied in detail since the mid 80s by several authors (Burrows & Lattimer bl86 (1986); Keil & Janka kj95 (1986); Pons et al. pns99 (1999, 2000, 2001)). The initial hot, lepton rich remnant, left behind following a successful core collapse supernova is known to evolve in a timescale of tens of seconds to form a cold ($`T<10^{10}`$ K) catalyzed neutron star. This process is usually classified in three main stages: i) the mantle contraction during which fast cooling of the outer regions takes place in about 0.5 s, with probably significant accretion; ii) deleptonization and consequently heating of the internal core as energetic neutrinos diffuse out leaving most of their chemical energy (analogously to Joule effect) on the way out; and iii) cooling by means of diffusion of (mostly) thermal neutrinos, resulting in a decrease of temperature from about 40-50 MeV to below 2-4 MeV, point at which the star becomes transparent to neutrinos.

Less is known about the early evolution of PNSs that are rapidly rotating, because a fully consistent study is a formidable task. It requires solving the neutrino transport equations in at least two-dimensions, coupled to general relativistic hydrodynamics that describes the fluid motion in rotating relativistic stars. In addition, the microphysical inputs also need to be treated carefully. As has been shown previously (Pons et al. pns99 (1999)), internal consistency between the equation of state of dense matter and the neutrino opacities describing the interaction between neutrinos and matter is needed to achieve reliable results.

A few attempts to treat the problem in a simplified way exist (Romero et al. rom92 (1992); Goussard et al. G97 (1997, 1998); Sumiyoshi et al. sumi99 (1999); Strobel et al. SSW99 (1999); Yuan & Heyl YH03 (2003)) based on the analysis of temporal sequences of quasi-equilibrium models of PNS, imposing some ad-hoc thermodynamical profiles (i.e., constant temperature, constant entropy and/or constant neutrino fraction). The assumption of quasi-equilibrium is well justified, since the hydrodynamical timescale ($`10^3`$ s) is much smaller than the timescale in which substantial thermodynamical changes occur (diffusion timescale $`1`$ s). However, it is unclear whether or not an isothermal or isentropic profile can effectively mimic reality. In addition, there are some physical constraints such as conservation of energy (the total gravitational mass must decrease consistently with the neutrino luminosity), or conservation of angular momentum, that need to be satisfied when one studies a temporal sequence of quasi-equilibrium models of PNS. These constraints help to reduce the parameter space as one calculates the thermodynamical evolution, although there is much uncertainty about the transport of angular momentum. In general angular momentum losses by neutrinos are supposed to be small (Kazana Kaz77 (1977)). Neutrinos could redistribute or take away a fraction of the initial angular momentum, but quantifying this amount needs of multidimensional transport simulations. It is equally controversial how important is turbulent transport, or the role of gravitational radiation. Because of these many unresolved issues it is usually assumed that the angular momentum is approximately conserved during the Kelvin-Helmholtz stage. In addition, the presence of large magnetic field could be crucial. It has been recently shown that turbulent mean–field dynamo action can be effective for PNSs with periods shorter than 1 s, which would generate very strong magnetic fields in the interior (Bonanno et al. BRU (2003)).

The whole problem is far to be solved, and will require some serious computational effort but, in the meantime, the intention of this paper is to begin to understand qualitatively what kind of outcome will be revealed by such numerical simulations, prior to engaging large scale simulations. In this line, we improve in several respects the few existing previous works. First, instead of isentropic or isothermal models, we use realistic profiles coming from 1D simulations, conveniently rescaling the temperature and chemical profiles as functions of density; second, we do not restrict ourselves to the rigid rotation case, as most of previous studies, and include differential rotation in our analysis; last but not least, we check that fixing the baryonic mass and the angular momentum results in some reasonable time dependence to the integrated relevant physical quantities such as the gravitational mass. In this way, we can calculate a complete evolutionary sequence of a PNS .

The plan of the paper is as follows. In 2 we briefly review the relativistic description of rotating stars, and the numerical method used to obtain equilibrium configurations. A similar brief description of the thermodynamics and the equation of state is presented in 3, since for both issues the details have been published elsewhere. Section 4 discusses the choice of parameters for the initial models. Our results are presented in section 5. In the first part we describe rigidly rotating PNSs, while the description of differentially rotating PNSs is given in 5.2. In both cases, we give all the relevant parameters for maximally rotating configurations as well as the evolution of more realistic sequences. Finally, our main conclusions and findings are summarized and discussed in 6.

## 2 Numerical code for stationary axisymmetric spacetime

Since the time scale for the matter and the gravitational field to settle in stationary equilibrium is much shorter than all other involved times111unless a gravitational dynamical instability occurs., a simplifying but plausible approximation consists of mimicking the time evolution of the spacetime by obtaining a collection of stationary, axisymmetric and asymptotically flat snapshots, all of them being solutions of the Einstein equations, but with a source term that depends on the snapshot.

This approximation has already been adopted by several authors in the past (Goussard et al. G97 (1997, 1998); Strobel et al. SSW99 (1999); Yuan & Heyl YH03 (2003)), our code being moreover very similar to the one of Goussard et al. (G97 (1997, 1998)), since we also employ a fully relativistic spectral code based on the work of Bonazzola et al. (BGSM (1993)). In the following, this article will be referred to as BGSM, and we address the interested reader to it for more details about the algorithms. Let us just remind that the BGSM approach uses the $`3+1`$ formulation of general relativity, maximal slicing-quasi-isotropic coordinates, with the assumption that matter is free of convective motions. The latter hypothesis is probably too strong, since PNSs are known to be convectively unstable (Keil et al. KJM96 (1996); Miralles et al. MPU (2002)), but it implies the circularity of the stationary axisymmetric spacetime (Carter c69 (1969)) and greatly simplifies the calculations.

Under these approximations, the metric is written 222For all this article, we are using natural units in which the speed of light is equal to unity.

$`\mathrm{d}s^2`$ $`\widehat{=}`$ $`g_{\mu \nu }\mathrm{d}x^\mu \mathrm{d}x^\nu `$ (1)
$`\widehat{=}`$ $`\left(N^2N_\phi N^\phi \right)\mathrm{dt}^22N_\phi \mathrm{dt}\mathrm{d}\phi `$
$`+`$ $`{\displaystyle \frac{A^4}{B^2}}\left(\mathrm{dr}^2+r^2\mathrm{d}\vartheta ^2+B^4r^2\mathrm{sin}^2\vartheta \mathrm{d}\phi ^2\right),`$
where we have introduced the usual notation of the $`3+1`$ formalism: $`N`$ is the lapse, $`N^\phi `$ is the third component of the shift $`3`$-vector, with $`N_\phi \widehat{=}g_{\phi i}N^i=A^4B^2N^\phi r^2\mathrm{sin}^2\vartheta `$ being the covariant $`\phi `$-component of the latter. Thanks to the hypothesis of axisymmetry and stationarity, the functions $`N,N^\phi ,A`$ and $`B`$ depend only on $`r`$ and $`\vartheta `$ .

Despite the general formalism and structure of the code is very similar to that of Goussard et al. (G97 (1997, 1998)), there are a few differences that should be mentioned:

* we use a more recent version of the spectral BGSM code (see Gourgoulhon et al. GHLPBM (1999)), in which at least two shells map the interior of the star, in order to better describe the non-spherical surface. This code is based on the C++ library LORENE333 http://www.lorene.obspm.fr, a software package for numerical relativity freely available under GNU license;
* the equation of state (EOS) is written in tabular form and interpolated using cubic Hermite polynomials, as defined in Nozawa et al. (NSGE (1998)), according to a general technique presented in Swesty (S96 (1996)) that preserves thermodynamical consistency (see also further remarks below);
* during the evolution of PNSs we do not restrict ourselves to 3 or 4 typical evolutionary snapshots that are just isentropic or isothermal profiles. Instead we adopt more realistic profiles, rescaled from the results of 1D simulations with neutrino transport (Pons et al. pns99 (1999), more details are given in the next section).

For clarity, let us be more specific about the last remark. From the hydrodynamical point of view, one must deal with the conditions for the existence of stationary solutions of the Einstein equations in presence of a rotating perfect fluid (see also BGSM and Goussard et al.). The main issue is that the relativistic Euler equations, derived from the conservation of the energy-momentum tensor projected into a space-like $`3`$-surface, are reduced in the case of stationary rotating motion to

$$\frac{_iP}{e+P}+_i\mathrm{ln}\left[\frac{N}{\mathrm{\Gamma }}\right]=F_i\mathrm{\Omega },$$
(2)
where

* $`P`$ is the pressure and $`e`$ the total energy density, both measured in the comoving frame;
* $`N`$ is the already defined lapse function, $`\mathrm{\Gamma }`$ the Lorentz factor linking the Eulerian observer and the comoving observer, $`\mathrm{\Omega }\widehat{=}u^\phi /u^t`$ the angular velocity, and $`F\widehat{=}u_\phi u^t`$ with $`u`$ being the $`4`$-velocity of the fluid;
* the index $`i`$ is either $`1`$ or $`2`$, since all the previous quantities depend only on $`r`$ and $`\vartheta `$.

As explained in BGSM Eq. (2) does not admit a first integral for arbitrary thermodynamics or rotation profiles $`\mathrm{\Omega }[r,\vartheta ]`$, but it does in some simple cases.

Concerning rotation, the right hand side of equation (2) is a first integral if either $`\mathrm{\Omega }`$ is constant in space (rigid rotation), or the function $`F`$ can be locally written as $`F=F[\mathrm{\Omega }]`$. The latter leads to the relation

$$F[\mathrm{\Omega }]\frac{A^4B^2r^2\mathrm{sin}^2\vartheta (\mathrm{\Omega }N^\phi )}{N^2A^4B^2r^2\mathrm{sin}^2\vartheta (\mathrm{\Omega }N^\phi )^2}=\mathrm{\hspace{0.17em}0},$$
(3)
between $`\mathrm{\Omega }`$ and the coefficients of the metric, which allows the determination of the profile of angular velocity.

For the differential rotation, we adopt the law proposed by Komatsu et al. (keh89 (1989)) and used by many authors (Komatsu et al. keh89 (1989); Goussard et al. G97 (1997, 1998); Baumgarte et al. Bau00 (2000))

$$F[\mathrm{\Omega }]=R_0^2(\mathrm{\Omega }_c\mathrm{\Omega }),$$
(4)
where $`\mathrm{\Omega }_c`$ is the limit of the function $`\mathrm{\Omega }[r,\vartheta ]`$ on the rotation axis and $`R_0^2`$ is a parameter with the dimension of the square of a length, noted $`A`$ in Komatsu et al. (keh89 (1989)). We shall discuss later our choices for the values of these parameters. It must be remarked that despite its simplicity, this law gives reasonable rotation profiles qualitatively similar to those obtained from core collapse simulations (see section 4).

Regarding the EOS, it can be shown that isentropic or isothermal profiles are sufficient conditions for Eq. (2) to be integrable. Yet another possibility which is less restrictive exists : assuming that the total energy density $`e`$ can be written as an effective function of the pressure $`e=e[P]`$. Indeed, in this case, even if the actual EOS is not barotropic, assuming that the EOS is a one–parameter function ensures the existence of the first integral

$$\widehat{H}+\mathrm{ln}\left[\frac{N}{\mathrm{\Gamma }}\right]+F[\mathrm{\Omega }]d\mathrm{\Omega }=\text{const.},$$
(5)
where the function $`\widehat{H}`$ is defined by

$$\widehat{H}\widehat{=}\frac{\mathrm{d}P}{e[P]+P},$$
(6)
with suitable boundary conditions. In this work, we have rescaled the temperature and composition profiles assuming, at each evolutionary time, that they depend only on the baryon density $`n_B`$. In this way, the EOS behaves effectively as a barotrope, and the pseudo-enthalpy function (6) can be written as

$$\widehat{H}\widehat{=}_0^{n_B}\frac{\mathrm{d}P[n_B]}{\mathrm{d}n_B}\frac{1}{e[n_B]+P[n_B]}dn_B,$$
(7)
with the effective functions $`P[n_B]`$ and $`e[n_B]`$ depending only on the baryonic density $`n_B`$.

At zero temperature, $`\widehat{H}`$ is equal to the logarithm of the specific enthalpy

$$H=\mathrm{ln}\left[\frac{e+P}{n_Bm_B}\right],$$
(8)
where $`m_B`$ is the baryon mass. As a basic test of the numerical method of integration, we have checked, in the case of tabulated EOS at $`T=\mathrm{\hspace{0.17em}0}`$, the agreement between $`\widehat{H}`$ calculated numerically and the analytical result (8). Due to the finite number of points in the tabular EOS the agreement is not complete, but the discrepancy was always found to be smaller than $`1\%`$ in all tested cases (even with strong differential rotation).

## 3 Equation of state and thermodynamics

Most of previous studies of structure of rotating neutron stars have only been performed with polytropes, or at most zero temperature EOS. This is not enough for the purpose of studying PNS evolution, since thermal and chemical effects are relevant and need a careful treatment. Only in a few works (Goussard et al. G97 (1997, 1998); Strobel et al. SSW99 (1999)) a realistic, finite temperature EOS has been used. The works by Goussard et al. (G97 (1997, 1998)) used a Skyrme-like potential model (Lattimer & Swesty LS91 (1991)), while the more recent work by Strobel et al. (SSW99 (1999)) used a Thomas-Fermi model with momentum and density dependent effective nucleon-nucleon interaction developed by Myers & Swiatecky (MS90 (1990)).

In this paper, we employ an EOS in the framework of relativistic field theory, in which the interaction between nucleons is mediated by the exchange of three meson fields. The coupling constants are chosen to reproduce some bulk properties of nuclear matter such as the nuclear saturation density, symmetry energy or nuclear compressibility. The formalism is widely used, and to avoid unnecessary repetition we refer to the complete review by Prakash et al. (Pra97 (1997)) and references therein. The particular choice of the coupling constants is detailed in Pons et al. (pns99 (1999)).

Figure 1 shows profiles from 1D simulations of a $`M_B=1.6M_{}`$ spherical star from Pons et al. (pns99 (1999)). In previous works (Goussard et al. G97 (1997, 1998); Strobel et al. SSW99 (1999)) the evolution of the star was approximated by choosing either simplified profiles of constant entropy or temperature, or the combination of a constant low entropy core ($`s=12`$) surrounded by a high entropy mantle ($`s=45`$) to mimic the situation in the early Kelvin-Helmholtz stage. It was shown that the choice of the critical density that separates the hot envelope from the inner core affected substantially the results. To reduce the uncertainties associated to the choice of this critical density and the fact that PNSs are not isothermal or isentropic configurations, we have chosen to use the thermodynamical structure from 1D simulations as a better approximation to reality, and assume that at each time the temperature and lepton fraction have the same dependence on the baryon density as the non-rotating models, $`T=T[n_B]`$ and $`Y_L=Y_L[n_B]`$. More explicitly: given the results from 1D simulations, the profiles of density, temperature, lepton fraction, entropy and pressure are tabulated at each time step, and these tables are taken as effective one–parameter EOSs where the independent variable is the baryon density ($`n_B`$) and the rest of variables are functions of only $`n_B`$.

Notice that, as we explained in the previous section, this (pseudo–barotropicity) assumption is a sufficient condition to ensure that the equation of motion has a first integral, making affordable the computation of the models. Of course, to know the exact profiles consistently, one should solve the neutrino transport equations in 2D along the whole evolution, but this requires a large technical and computational effort that is out of the scope of this paper.

## 4 Choice of the initial models

One of the first issues that must be addressed to study the early evolution of neutron stars is the initial model. In Pons et al. (pns99 (1999)) it was shown that for spherical models the overall evolution is not very sensitive to modifications of the initial thermodynamical profiles, while the total mass of the star is the parameter that mostly affects the subsequent evolution. It seems therefore natural to chose a canonical star with a baryon mass 1.6 $`M_{}`$, that corresponds to a gravitational mass of the old, cold configuration of about 1.44 $`M_{}`$. For completeness, we have also studied a model with a baryonic mass of 1.2 $`M_{}`$. Notice that one cannot use the effective barotropes obtained for a PNSs with a given mass in the calculations of sequences with a different mass, because the thermodynamic properties are different.

More difficult is to guess the particular rotation properties of a newly born neutron star. Newly born neutron stars are expected to rotate non rigidly, as opposite to old neutron stars in which the various viscous mechanisms had enough time to act. Even in the case that the iron core of the Supernova progenitor was rigidly rotating (a quite reasonable assumption, since it is an old massive star), the process of collapse would generate a significant amount of differential rotation. Indeed, recent calculations of stellar evolution indicate that the iron core of massive stars is rotating almost rigidly (Heger et al. Heg00 (2000, 2003)) and this is the approximation used in the most recent simulations of rotating stellar core collapse (Müller et al. Ewa03 (2003)). Based on this works, we have decided to study a few models from Dimmelmeier et al. (DFM (2002), DFM hereafter) and to run some more simulations of similar models of stellar core collapse varying the total initial angular momentum. Our models correspond to those labeled by A1 in DFM, i.e., the iron core of the progenitor is almost rigidly rotating. In Figure 2, we show profiles of the angular velocity (upper panel) and density (lower panel) of the inner 150 km, as a function of the equatorial radius, immediately after the PNS is formed. The different lines correspond to models with a different amount of initial angular momentum, with values of the ratio between rotational to gravitational potential energy ($`|T/W|`$) of 0.9% (models B3 in DFM), 0.5% (B2), 0.25% (B1) and 0.05%, respectively. We included the last lower value because it is more realistic (Heger et al. Heg03 (2003); Müller et al. Ewa03 (2003)) than those studied in DFM. In the figure we also show with dots the simple fit

$$\mathrm{\Omega }[r,\vartheta ]=\mathrm{\Omega }[r\mathrm{sin}\vartheta ]=\frac{\mathrm{\Omega }_cR_0^2}{R_0^2+r^2\mathrm{sin}^2\theta }$$
(9)
which corresponds to the Newtonian limit of $`\mathrm{\Omega }`$ consistent with Eq. (4). In this limiting case, we now see that $`R_0`$ is the value of the radius at which, in the equatorial plane ($`\vartheta \widehat{=}\pi /2`$), the angular velocity is half its limiting central value, $`\mathrm{\Omega }_c`$. Even if this remark is valid only in the Newtonian case, it is nicely verified for relativistic stars as well. For all four models we have taken $`R_0^2=50`$ km2, while $`\mathrm{\Omega }_c`$ takes the values of 4500, 3500, 2600, and 1300 rad/s, from top to bottom. The agreement between the simple law (9) and the results from simulations is quite acceptable. In all the models the spatial scale of variation of the angular velocity ($`R_0`$) is the same because the density profiles were very similar (see Figure 2, bottom). Varying the EOS might lead to a different density distribution, in such a way that one should expect smaller $`R_0`$ for softer EOSs (more compact PNS), and larger $`R_0`$ for stiffer EOSs. We must remark again that our assumption of having stationary axisymmetric solutions implied that the function $`F`$ on the r.h.s. of the equation of stationary motion (2) is only a function of the angular velocity. In the Newtonian case, this is equivalent to say that the angular velocity depends only on the distance to the axis ($`r\mathrm{sin}\vartheta `$). In the relativistic case, there are small corrections but the overall distribution of angular velocity is nearly cylindrical, as can be seen in Fig. 3.

In reality, immediately after core collapse and bounce, there is no reason to expect that the distribution of angular velocity corresponds to that of stationary equilibrium, and it may take several rotation periods for the star to relax to some stationary solution. A more detailed analysis of all these issues is out of the scope of this paper, and requires a parametric study of rotating core collapse simulations. However, for our purpose, which is just to understand in which range of parameters one should move when using Eq. (4), the conclusion is that a reasonable realistic model with differential rotation should be consistent with $`R_0`$ of the order of 10 km. Considering that there might be a number of processes that could lead to rigid rotation on a dynamical timescale, such as shear or magneto-rotational instabilities, we will study both cases, with substantial differential rotation and with rigid rotation. Notice that previous works (Goussard et al. G98 (1998)) on differentially rotating PNSs have focused on the parameter space region corresponding to $`R_01`$ km, which seems too low according to our results.

A final remark concerns the geometry of the PNSs. Although for progenitors with moderate angular momentum, conservation of angular momentum leads to a characteristic profile of angular velocity that decreases with increasing distance to the rotation axis, in the case of very fast rotation, or unrealistic large differential rotation, the usual oblate shape of the PNS can become toroidal, with a maximum in density off axis. In this case the profile of the angular velocity can exhibit also a maximum at some distance from the axis. We will discuss below for which models a toroidal shape is found, but unfortunately our numerical algorithm does not allow us to compute equilibrium configurations with zero density at the origin.

## 5 Results

We have considered a number of thermodynamical profiles corresponding to snapshots of the early evolution of PNSs, and we have calculated the rotational properties in two different cases: rigidly rotating and differentially rotating PNSs. The left part of Table 1 shows some properties (central baryonic density, gravitational mass, and radius) of the non rotating models with a baryonic mass of 1.6 $`M_{}`$. The models are labeled according to the evolutionary time (in s after birth) in the non-rotating case. Notice that for rotating PNSs the evolution time labeling the models is only an indication, since in a self-consistent simulation one would expect faster evolution. The reason of this faster evolution is that the lower densities, and consistently lower temperatures because the entropy is conserved, result in larger neutrino mean free paths. Therefore the neutrino transport would be more effective.

In order to investigate the influence of the low density ($`<10^{11}`$ g/cm3) EOS, we have also substituted the zero temperature BPS equation of state (Baym et al. bps (1971)) by the LS (Lattimer & Swesty LS91 (1991)) EOS in some models at early times. The results for the models with the LS EOS are shown between parenthesis in the first four rows of Table 1. At later times, after the mantle loses most of its thermal and lepton content and contracts, the influence of the low density EOS becomes less important and both EOSs give essentially the same result, because the amount of mass and the thickness of the low density layer is rather small. The left hand side of Table 2 shows the equivalent information but for a non rotating NS with a total baryonic mass of 1.2 $`M_{}`$.

### 5.1 Rigidly rotating PNSs

In the right hand side of Tables 1 and 2 we summarize the properties of maximally rotating PNSs following the same thermodynamical sequence as in the non-rotating case. Focusing at early times ($`t<0.5`$ s) we see that the central densities of non-rotating and maximally rotating PNSs are not terribly different. The reason is that these objects are only slightly bounded and they cannot rotate very fast, reaching the mass shedding limit at angular velocities of hundreds of Hz. Therefore, the rotational energy is not large and the models are not very different from the non-rotating case. This can also be seen in the values of the quantity $`|T/W|`$, which is about 0.03–0.05 for the BPS EOS and 0.06–0.07 when using the LS EOS at low density. The differences between these two EOSs come from the fact that the BPS is at zero temperature while the LS is at finite temperature. Due to the more extended mantle, the total angular momentum is systematically large, of about 30$`\%`$, for the models with finite temperature at $`\rho <10^{11}`$ g/cm3.

After about 0.5-1.0 seconds, the mantle has contracted and lost its thermal and lepton content, so that the PNS becomes more compact and bound. This can be deduced by looking at the central density, which exhibits a sudden increase at the end of this first epoch. For the non-rotating models it happens between 0.5 and 1.0 seconds, while for the maximally rotating models it seems to be delayed until some time between 1.0 and 2.0 seconds, due to the role of the centrifugal force. For these more compact configurations, the limiting angular velocities are much larger (3000-6000 rad/s), and their central densities are considerably lower than those of the non-rotating PNSs. The ratio between polar and equatorial radii of maximally rotating PNSs turns out to be rather insensitive to the thermodynamics, being about 0.65 for early models and decreasing to 0.55 for late configurations. With respect to the $`|T/W|`$ ratio, its maximum value is about 0.1 at the end of the evolution, still below the threshold to undergo dynamical instabilities. As a rule of thumb, we found that the maximum angular velocity can be estimated with good accuracy with the simple fit $`\mathrm{\Omega }_K0.58\sqrt{GM/R^3}`$, valid for both masses (1.6 and 1.2 $`M_{}`$) and all thermodynamical profiles studied. Despite our models are not close to the maximum mass for a given EOS, this result is in agreement with the empirical formula derived by Haensel & Zdunik (HZ89 (1989)) relating the maximum rotation frequency to the mass and radius of the maximum mass static configuration. It seems that this empirical relation is valid in general.

The evolution of the angular momentum deserves a separate discussion. The first result is that including thermal effects at low density (LS) allows for significantly larger total angular momentum. However, this discrepancy disappears rapidly, and the maximum angular momentum decreases quite fast and monotonously with time, which has evolutionary implications. Whether or not neutrino transport can redistribute and take away a substantial fraction of the initial angular momentum is unclear and needs of multidimensional transport simulations to be clarified. In general, it is assumed that the angular momentum is approximately constant during the Kelvin-Helmholtz epoch, and it only varies on timescales larger than the diffusion timescale. If this was the case, PNSs born with maximal angular momentum cannot radiate it via neutrino diffusion with the same efficiency as the binding energy, and this may result in significant mass loss after reaching the mass shedding limit. Other possibility is that PNSs are not born maximally rotating, the angular momentum is (almost) conserved in a timescale of seconds, and the star speeds up as it contracts and evolves thermodynamically. We have simulated this scenario by taking a sequence of configurations with constant angular momentum ($`J=1.5GM_{}^2/c`$), always with fixed baryonic mass. Results are shown in Figure 4. It can be seen how a model that initially had a $`|T/W|`$ ratio of only 0.02 speeds up in a timescale of a second and then slowly approaches the mass shedding limit in a few seconds. We expect that fast (but not maximally) rotating PNSs might reach the threshold of some types of instabilities (CFS, bar–mode, …) on a diffusion timescale, say within 10 s after birth. It will depend on the efficiency of neutrinos to carry away angular momentum, which is currently under investigation.

The above discussion was focused on results for the canonical NS with a baryonic mass of 1.6 $`M_{}`$. For completeness we also give results for a model with a lower baryonic mass of 1.2 $`M_{}`$ in Table 2. The results are qualitatively similar to those of Table 1, with only some quantitative differences. Accordingly to the lower mass, the maximum angular velocities and $`|T/W|`$ ratio are also lower.

### 5.2 Differentially rotating PNSs

As discussed in the previous section we will restrict to realistic values of $`R_0`$ in the range from 10 km (significant differential rotation) to 50 km (almost rigid rotation). A different scenario would be a hot NS resulting from the merger of two neutron stars. In this case, the distribution of angular momentum is more complex and depends not only on the initial total angular momentum, but also on a number of different parameters as the mass ratio or the relative orientation of the spins of the two stars. Lower values of $`R_0`$ might be reached in this latter scenario, but above all, the rotation profile would probably be very different from Eq. (4), which is in agreement with a stationary configuration, the case of a merging being far from a stationary situation.

As before, the baryonic mass of the star is conserved and fixed to 1.6 $`M_{}`$ for all models. Within the range of parameters mentioned above, we have calculated a large number of models of differentially rotating PNSs with either a fixed $`R_0`$ of 10, 20, 50, and 100 km or varying $`R_0`$ according to the circumferential equatorial radius of the star ($`R_0=\mathrm{\hspace{0.17em}0.5}R_{eq}`$). In Figures 5 and 6 we show iso-enthalpy surfaces for two representative models at an evolutionary time of $`t=1`$ s. Despite both configurations have similar equatorial radii ($`30`$ km), they are qualitatively very different. In Figure 5 the parameter $`R_0=100`$ km, i.e., the model has a very low degree of differential rotation. Even very close to the mass shedding limit (0.995 of the maximum angular velocity) the star shows a regular oblate profile. On the contrary, the model in Figure 6 corresponds to the case of significant differential rotation ($`R_0=10`$ km). With strong differential rotation we have found that above some critical angular velocity (950 rad/s, for this model) the shape of the star becomes toroidal, with a maximum density off axis.

The most relevant results are presented in a compact form in Figures 7 and 8, where we show the total angular momentum and the rotation parameter $`|T/W|`$ as a function of the central angular velocity $`\mathrm{\Omega }_c`$, for different values of $`R_0`$. For comparison, we have also included the results for rigid rotation (solid curve), together with those for $`R_0=50`$ km (dash-dotted line), $`R_0=20`$ km (thick dashes), $`R_0=R_{eq}/2`$ (thin dashes) and $`R_0=10`$ km (dots). For $`R_0100`$ km, the results are so close to the rigid rotation case that we do not show them in the figures.

The thermodynamical profiles for each panel are those corresponding to evolutionary times of the non-rotating models of $`t=0.5`$ (top), $`1.0`$ (middle) and 10 s (bottom panel) after formation. From the results of the figure we can see that for $`R_0`$ of the order or larger than 50 km the effects of differential rotation are unimportant. Only for $`R_0<20`$ km and $`R_0=0.5R_{eq}`$ the differences are significant. As discussed in 4, results from axisymmetric core collapse seem to indicate that the scale in which one should expect variations of the angular velocity is of about 10 km, therefore the last two cases may be looked at as the most realistic, with the caveat of being constrained to the particular rotation law we used.

In the model corresponding to $`t=1`$ s it is visible a sudden increase of the angular momentum at a critical value of the angular velocity, different for each value of $`R_0`$. This is an effect of the phase transition from bulk matter to nuclei, and it is associated to the point when the angular velocity becomes so large that the central density is below nuclear saturation density ($`\rho _0`$). This transition is not observed at earlier times because matter at low density is still very hot, and the thermal effects suppress or smooth out the kink. At later times, the star is more compact and the central density is always above $`\rho _0`$, even close to the mass shedding limit. For the models with strong differential rotation, we had to stop our calculations before reaching the mass shedding limit because of numerical limitations associated to the special toroidal geometry that arises (see Fig. 6). These cases are marked by an asterisk in Fig. 7 and 8. In principle, there could be configurations with slightly larger angular momentum.

An expected result is that the central angular velocity can be up to a factor 5 to 10 larger in the case of differential rotation, but the maximum angular momentum is only about a fifty percent larger than for rigid rotation. More important is the variation in the maximum of $`|T/W|`$. While for rigidly rotating stars it is about 0.1, for PNSs with a large degree of differential rotation it can be as high as 0.2. The fact that differential rotation allows for the existence of equilibrium models at large $`|T/W|`$ have also been found by Baumgarte et al. (Bau00 (2000)) and Yoshida et al. (YRKE (2002)) for polytropes in the relativistic Cowling approximation using the same rotation law. A maximum value of $`|T/W|=0.2`$ is not enough to reach the dynamical instability threshold ($`0.27`$), but it is sufficient to allow either the secular, gravitational wave driven, instability (Chandrasekhar cha (1970); Friedman & Schutz FS (1978)) that happens at $`0.10.15`$ or the recently proposed low $`|T/W|`$ dynamical instabilities (Shibata et al. Shi02 (2002, 2003); Watts et al. WAJ (2003)). Notice that despite our calculations cannot be carried up to the mass shedding limit for some models, we are not far from it and the numbers are close to the maximum value. This can be seen in Figure 9, when we show the polar to equatorial axis ratio as a function of $`\mathrm{\Omega }_c`$. For rigid rotation, or moderate differential rotation the calculations end at the Keplerian limit, while for strong differential rotation the calculation ends because of numerical problems (marked with asterisks). The slope of the curves is becoming more and more vertical, and the axis ratio is as low as 0.4, smaller than the minimum in the rigid case (0.6). This is an indication that the shape of the star is changing fast and only slightly larger angular velocities can be reached before the star starts loosing mass.

The next question one can pose is whether or not dynamical instabilities arise in a real case. As for rigidly rotating stars, we now impose that the angular momentum is fixed along the evolution, as well as the baryonic mass, and study how the different variables evolve under such constraints. Figure 10 shows the temporal evolution of the central angular velocity (top) and of the $`|T/W|`$ ratio (bottom) of a PNS differentially rotating with $`R_0=10`$ km and a fixed angular momentum of $`J=1.5GM_{}^2/c`$. As in the case of rigid rotation, $`|T/W|`$ increases as the star contracts and looses its binding energy in a neutrino diffusion timescale of about 10 s. Thus, if there are no significant losses of angular momentum, the various types of dynamical instabilities associated to critical values of the rotation parameter may arise not at the very beginning, but several seconds after formation.

As a consistency check of our assumptions and many simplifications, one can make a simple estimate of the neutrino luminosity considering that the variation in the gravitational mass of the PNS is considered to be equal to the luminosity,i.e., $`L_\nu =\frac{dM_G}{dt}`$. The result is shown in Figure 11. Both the order of magnitude and the exponential decay of the luminosity in a timescale of 10 seconds are consistent with the results of detailed simulations with neutrino transport in spherical symmetry (Pons et al. pns99 (1999)). It is interesting that the difference between the two models studied with $`R_0=10`$ km (solid) and $`R_0=20`$ km (dashes) is only visible after 10 seconds of evolution, when the stars become more compact. Therefore, imposing conservation of angular momentum and baryonic mass, and building a sequence keeping fixed $`R_0`$ is not violating any basic physical law, such as conservation of energy, and can be taken as a qualitative approach to the real case.

## 6 Conclusions

We have approached the problem of the evolution of rotating protoneutron stars by constructing evolutionary sequences of axisymmetric stationary configurations in General Relativity. The thermodynamical structure and evolution have been taken and extrapolated from simulations in spherical symmetry that included neutrino transport. Although this is a crude simplification, it already gives an interesting insight about how the different relevant quantities can evolve as the rotating PNS loses its lepton content and its excess binding energy, and contracts. Moreover, we have found that the luminosity estimates are not terribly different from what one expects.

A special effort has been made to understand in which space of parameters we should be in a realistic case. The biggest uncertainty concerns the rotation law that PNS have at birth. By analyzing results from simplified core collapse simulations, it seems that a typical scale for variations of the angular velocity is about 10 km, and that conservation of angular momentum during the collapse of a stellar core (initially rigidly rotating) does not seem to allow for angular velocities varying in a length–scale shorter than a few km. Less is known about the angular distribution, except that the most recent simulations show the presence of important meridional currents and some turbulent motion. For simplicity, we restricted ourselves to the stationary case. Stationarity implies a quasi-cylindrical distribution (with deviations due to relativistic corrections) of the angular velocity. This stage only can be reached after several dynamical and rotation periods, after the PNS had time to relax. Therefore one must be aware that the first 0.5 s are probably far from stationarity, but after that evolution proceeds in a quasi-stationary way, except for low velocity convective motions. From our study of quasi-stationary sequences we can draw a few interesting qualitative results.

i) For rigidly rotating stars, the mass shedding limit is well approximated by the simple law $`\mathrm{\Omega }_K0.58\sqrt{GM/R^3}`$, despite complications in the EOSs. This empirical formula is valid for stars of different masses within a 5% accuracy for the EOSs we considered. We believe that this result is general and would be valid for other EOSs or stars with different masses.

ii) For each instant in the evolution, stars with strong differential rotation can have 5 to 10 times larger central angular velocities, and accommodate about a fifty percent more angular momentum. The maximum specific angular momentum $`J/M`$ varies between $`(12)GM_{}/c(0.51)\times 10^{16}`$ cm2/s, depending on the degree of differential rotation.

iii) When $`R_020`$ km there is a critical value of the angular velocity above which the shape of the star becomes toroidal, with a maximum density off axis. We found that this configurations are stationary solutions of the hydrodynamics equations, but it is unclear under which circumstances this geometry is stable against perturbations. This issue deserves a separate study.

iv) The maximum value of $`|T/W|`$ obtained in the case of differential rotation is about 0.2, while for rigid rotation this is $`0.11`$. Thus differentially rotating PNSs might undergo the CFS instability, which arises at $`0.14`$ and, in any case, the recently discussed low $`|T/W|`$ instability (Shibata et al. Shi02 (2002, 2003); Watts et al. WAJ (2003)) is plausible to happen.

v) More interestingly, we found several situations in which, even if the initial model is not close enough to the critical value of $`|T/W|`$, as the star contracts in a neutrino diffusion timescale of 5-10 s, it speeds up entering the window of instability. An observational evidence of this effect could be a temporal shift of a few seconds between the neutrino luminosity peak and a gravitational wave burst in the event of a galactic Supernova. Ultimately, this depends on the initial amount of angular momentum, which is approximately equal to the angular momentum of the iron core of the progenitor. Recent stellar evolution calculations suggest that the specific angular momentum of the inner 1.7 $`M_{}`$ of a 15 $`M_{}`$ star can be as high as $`3\times 10^{16}`$ cm2/s if magnetic braking is neglected, or $`10^{15}`$ cm2/s if magnetic torques are included in the evolution (Heger et al. Heg03 (2003)). This corresponds to $`J/M`$ in the range $`J/M=(0.26)GM_{}/c`$. If the angular momentum happens to be in the upper region ($`J/M>2`$), centrifugal forces would stop the collapse before the PNS is formed. Intermediate values ($`J/M=1`$) may result in the formation of a rapidly rotating PNS that enters the instability region several seconds after birth. If magnetic braking is very effective during the evolution of a massive star, $`J/M<0.5`$ and the PNS will be formed after collapse without reaching extreme values of the angular velocities and $`|T/W|`$.

The next natural step to improve this work is to include the possible mechanisms to transport angular momentum between the different layers of the star, that may involve neutrino transport, turbulent transport, magnetic fields, neutrino viscosity, convective motion and/or angular momentum losses by gravitational wave emission. Unless the star is born with almost maximal angular velocities, since the different instabilities can arise in a timescale of seconds some of all of these dissipative mechanisms can modify our current vision of PNS evolution.

## Acknowledgments

We are grateful to H. Dimmelmeier for providing us with rotation profiles that we have used to constrain the parameters of the rotation law. We thank L. Rezzolla, J. Novak, S. Yoshida and D. Gondek-Rosinska for useful discussions and comments. This work has been supported by the EU Programme ‘Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137) and the Spanish Ministerio de Ciencia y Tecnologia grant AYA 2001-3490-C02. JAP is supported by a Ramón y Cajal contract from the Spanish MCyT and LV benefited from the Jumelage PAN-CNRS Astronomy France-Poland.


text:
# Pulsating White Dwarfs in Cataclysmic Variables: The Marriage of ZZ Cet and Dwarf Nova

## 1. Introduction

The DA white dwarf non-radially pulsators stars, known as ZZ Cet stars, have traditionally been found as isolated field objects. Until recently 34 were known (Fontaine et al. 2003), but a further tranche of 34 were just announced (D. Sullivan, this conference). The ZZ Cet pulsation strip in the HR diagram roughly extends over the range 11 000 K $`<`$ $`T_{eff}`$ $`<`$ 12 000 K (Bergeron et al. 1995), so it is of interest to see whether cataclysmic variables (CVs) with temperatures in that range also show pulsations (though the width and centre of gravity of the instability strip may be modified in the presence of accretion: Townsley & Bildsten 2003).

Until recently only one such hybrid CV/ZZ system was known: GW Lib (Warner & van Zyl 1998), which has a measured $`T_{eff}`$ of 14 700 K (Szkody et al. 2002a) – though this may be affected by non-allowance for the mysterious emission source that adds a 2.08 h photometric modulation independent of the 1.28 h orbital period (Woudt & Warner 2002). The principal periods in GW Lib are near 650 s, 376 s and 236 s.

A second CV/ZZ system, SDSS J161033.64-010223.3 (SDSS 1610 hereafter), was discovered in June 2003 (Woudt & Warner 2003) and was selected as a candidate on the basis of its spectrum as published in the first release of CVs in the Sloan Digital Sky Survey (Szkody et al. 2002b). SDSS 1610 resembles GW Lib in clearly showing absorption lines of the underlying white dwarf primary, as well as the emission lines characteristic of an accreting system. Its light curve is shown in Figure 1 (top panel). Its periodicities are near 607 s, 345 s and 221 s, with a harmonic at 304 s.

There is evidently a window of opportunity among the CVs of low rate of accretion ($`\dot{M}`$): the $`T_{eff}`$ of the white dwarf is determined by $`\dot{M}`$, largely through compressional heating in the interior (Sion 2003), and it happens that an $`\dot{M}`$ that maintains $`T_{eff}`$ in the instability strip is sufficient to produce Balmer emission lines but not too large to give the accretion disc a luminosity large enough to hide the flux from the white dwarf. From the observed depth of, e.g., H$`\beta `$, the fraction of flux contributed by the white dwarf can be estimated; for GW Lib this is $``$ 50% and in SDSS 1610 it is even greater.

The low $`\dot{M}`$ in these systems automatically leads to the expectation that they will be dwarf novae of very long outburst interval $`T_{out}`$ (see, e.g., Warner (1995)). Comparison with Z Cha, which has $`T_{out}`$ 50 d, shows that although the latter’s primary is visible in the spectrum, it is far more covered by accretion disc flux than in the known CV/ZZ stars. At the other end of the $`T_{out}`$ range, GW Lib has only had one observed outburst (in 1983) and SDSS 1610 has had none. Because the low $`\dot{M}`$ systems are intrinsically faint they will in general be apparently faint – it is the ability of the SDSS to find CVs to faint limits that has opened the possibility of increasing the number of known CV/ZZ stars.

## 2. Pulsation frequencies

Some of the ZZ Cet stars show frequency patterns that have been ascribed to direct resonance between the principal driving modes and the wealth of other available modes (O’Donoghue, Warner & Cropper 1992). Thus in VY Hor, GD 154 and PG 1351+489 the eigenfrequencies are almost completely described by the sequence $`nf`$ and $`(m+\frac{1}{2}+ϵ)f`$, where $`n`$ = 1, 2, … and $`m`$ = 0, 1, 2, …, and $`ϵ`$ is a small quantity (O’Donoghue, Warner & Cropper 1992; Robinson et al. 1978; Winget, Nather & Hill 1987). For VY Hor $`ϵ`$ = 0.037, for GD 154 $`ϵ`$ = –0.03 and for PG 1351 $`ϵ`$ = –0.03. In GW Lib and SDSS 1610, with the still rather limited observational data for the latter, we find a similar situation, but with the sequences $`nf`$ and $`(m+\frac{3}{4}+ϵ)f`$, where $`ϵ0.03`$ for GW Lib and $`ϵ0.00`$ for SDSS 1610. This may be a coincidence, but it could indicate that a different resonance condition is operating in accreting ZZ Cet stars.

## 3. Further discoveries

The publication of the second SDSS CV release (Szkody et al. 2003) revealed at least three further candidate CV/ZZ stars. We have observed SDSS J013132.39-090122.3 and SDSS J220553.98+115553.7 (see also Figure 1) and find that indeed both have ZZ Cet primaries. At the time of writing this leaves the strong candidate SDSS 1238 to be examined later in the observing season.

The Fourier transforms in Figure 2 show that in SDSS 0131 the dominant periods are near 595 s and 335 s with weak oscillations also near 260 s; SDSS 2205 has strong oscillations near 575 s and 330 s, with weak oscillations near 475 s. In both cases the ratio of the stronger periods is $``$ 1.75, as in GW Lib and SDSS 1610.

The total number of CVs in the first two SSDS releases is 60, of which we have found that probably $``$ 8% are CV/ZZ combinations. It is estimated that the final total of CVs found in the SDSS will be $``$ 400 (Szkody et al. 2002b), so in a few years we may expect to have more than 30 CV/ZZ stars to study in detail.

### Acknowledgments.

BW is supported by funds from the University of Cape Town; PAW is supported by strategic funds made available to BW by the University of Cape Town and by the National Research Foundation.

## References

Bergeron, P., Wesemael, F., Lamontagne, R., Fontaine, G., Saffer, R.A., & Allard, N.F. 1995, ApJ, 449, 258

Fontaine, G., Bergeron, P, Billères, M., & Charpinet, S. 2003, ApJ, 591, 1184

O’Donoghue, D., Warner, B., & Cropper, M. 1992, MNRAS, 258, 415

Robinson, E.L., Stover, R.J., Nather, R.E., & McGraw, J.T. 1978, ApJ, 220, 614

Sion, E.M. 2003, Proc. 13th European Workshop on White Dwarfs (Naples)

Szkody, P., Gaensicke, B.T., Howell, S.B., & Sion, E.M. 2002a, ApJ, 575, 79

Szkody, P., et al. 2002b, AJ, 123, 430

Szkody, P., et al. 2003, AJ, 126, 1499

Townsley, D.M., & Bildsten, L. 2003, ApJ, in press

Warner, B. 1995, Cataclysmic Variable Stars, Cambridge University Press

Warner, B., & van Zyl, L. 1998, IAU Symp. No. 185, 321

Winget, D.E., Nather, R.E., & Hill, J.A. 1987, ApJ, 316, 305

Woudt, P.A., & Warner, B. 2002, Ap&SS, 282, 433

Woudt, P.A., & Warner, B. 2003, MNRAS, in press


text:
# Pitfalls of Goodness-of-Fit from Likelihood

## I INTRODUCTION

> For every complex problem, there is a solution that is simple, neat, and wrong. H.L. Mencken

The complex problem considered here is goodness-of-fit (g.o.f.) for unbinned maximum likelihood fits in cases when binned g.o.f. methods and Kolmogorov-Smirnov are not well suited:

A physicist, having fit a complicated model to his multi dimensional data to obtain estimates of the values of certain parameters, is also expected to check how well the data match his model. In the sections that follow, we discuss a g.o.f. method, still occasionally used in high energy physics (HEP), that is simple, neat, and wrong.

## II THE SNW11footnotemark: 1 METHOD

22footnotetext: Simple, Neat, Wrong.

We start with a brief description of the method. (A true derivation, for obvious reasons, is not available.)

Maximum likelihood fits are performed by maximizing the likelihood $`L(\stackrel{}{\theta },\stackrel{}{x})`$ with respect to the (unknown) parameters $`\stackrel{}{\theta }`$ for fixed data $`\stackrel{}{x}`$.

Thus, the value of the likelihood provides the g.o.f. between the data and the probability density function (p.d.f.): The value of the likelihood at the maximum,

$$L_{\mathrm{max}}=L(\stackrel{}{\widehat{\theta }},\stackrel{}{x})$$

corresponds to the best fit—the smaller the likelihood, the worse the g.o.f., …

To calculate this “g.o.f.” P-value, we need the distribution of $`L_{\mathrm{max}}`$ for an ensemble of random $`\stackrel{}{x}`$ deviates from the p.d.f. using the true (but unknown) parameters $`\stackrel{}{\mathrm{\Theta }}`$.

We approximate this by replacing $`\stackrel{}{\mathrm{\Theta }}`$ with the parameter estimate obtained from the fit to the actual data.

This method has a long history of use in high energy physics. It’s recommended by several excellent statistical data analysis texts written by (and for) high energy particle physicists. Consequently, and because the method is “obvious”, it’s still being used in (some) HEP analyses.

Reference , written by a statistician and four physicists, describes the method, but criticizes:

> The likelihood of the data would appear to be a good \[g.o.f.\] candidate at first sight. Unfortunately, this carries little information as a test statistic, as we shall see…

Since this was ignored, maybe its warning was not strong enough. I have found no mention of the method in texts written (solely) by statisticians.

## III A SIMPLE TEST OF THE METHOD

> Always test your general reasoning against simple models. John S. Bell

Reference , following the above advice, tests the method against the p.d.f.

$$\frac{1}{\tau }e^{t/\tau }(t0)$$

where $`t`$ (we have in mind the decay-time of a particle) follows an exponential distribution, and $`\tau `$ (the mean lifetime) is a parameter whose value, being unknown, is estimated from data. The likelihood for $`N`$ observations $`t_i`$ is given by

$$\mathrm{ln}L=\underset{i=1}{\overset{N}{}}\left[\mathrm{ln}\tau +\frac{t_i}{\tau }\right]$$

The value ($`\widehat{\tau }`$) of $`\tau `$ that maximizes the likelihood, and the value ($`L_{\mathrm{max}}`$) of the likelihood at its maximum, are given by

$$\widehat{\tau }=\frac{1}{N}\underset{i=1}{\overset{N}{}}t_i\mathrm{ln}L_{\mathrm{max}}=N(1+\mathrm{ln}\widehat{\tau })$$

### III.1 The First Surprise

The value of the likelihood at its maximum (in this test case) is just a simple function of $`\widehat{\tau }`$—all samples with the same mean obtain the same “g.o.f.” value. This is a disaster for g.o.f. Even if the true value of $`\tau `$—call it $`𝒯`$—were known in advance, so that we could calculate the P-value associated with the observed $`\widehat{\tau }`$, merely comparing the $`\widehat{\tau }`$ of the data with $`𝒯`$ is not sufficient to show that the observed data are modeled well by the exponential distribution.

### III.2 The Second Surprise

Since under this method, our P-value ensemble is actually based on the value of $`\widehat{\tau }`$ computed from the data (not knowing the true value $`𝒯`$), we *always* obtain a P-value of about 50%, *for any data whatsoever*. This is a second disaster for g.o.f. By construction, the distribution of $`L_{\mathrm{max}}`$ from our ensemble of $`N`$-event pseudo experiments tracks the $`L_{\mathrm{max}}`$ observed from the data.

The fact that the method yields “reasonable” P-values has undoubtedly contributed to its longevity in practice: P-values very near 0 or 100% would have triggered further investigation.

### III.3 Lessons Learned

In this example, g.o.f. is equivalent to testing the single hypothesis: “The data are from an exponential distribution of unspecified mean.” $`L_{\mathrm{max}}`$ provided no information with respect to this hypothesis.

What went wrong? In our test case, the likelihood could be expressed as a function of just the parameter and its maximum likelihood estimator (m.l.e.): $`L(\tau ;\widehat{\tau })`$. *All* data samples with the same m.l.e. gave the same “g.o.f.”

Exactly the same thing happens in the Gaussian (normal) case—the likelihood can be written using solely the 2 parameters and their estimators: $`L(\mu ,\sigma ;\widehat{\mu },\widehat{\sigma })`$.

Other “textbook” distributions—scaled gamma, beta, log-normal, geometric—also fail in the same way. Geometric is a discrete distribution, so the problem is not restricted to the continuous case.

## IV MORE TROUBLE: NON INVARIANCE

Returning to our exponential example, suppose we make the substitution $`t=x^2`$. The p.d.f. transforms as

$$\frac{1}{\tau }e^{t/\tau }dt=\frac{2x}{\tau }e^{x^2/\tau }dx$$

and the g.o.f. statistic is now calculated as

$$\mathrm{ln}L=\underset{i=1}{\overset{N}{}}\left[\mathrm{ln}\tau +\frac{x_i^2}{\tau }\mathrm{ln}(2x_i)\right]\widehat{\tau }=\frac{1}{N}\underset{i=1}{\overset{N}{}}x_i^2$$

$`\mathrm{ln}L_{\mathrm{max}}`$ $`=`$ $`N(1+\mathrm{ln}\widehat{\tau }){\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{ln}(2x_i)`$
$`=`$ $`N(1+\mathrm{ln}\widehat{\tau }){\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{ln}(2\sqrt{t_i})`$
That is, the “g.o.f.” statistic is not invariant under change of variable in the continuous p.d.f. case. (The value of the m.l.e. is, of course, invariant.)

Under change of variable, the “g.o.f.” statistic picks up an extra term from the Jacobian—an extra function of the data. We’re free to choose any transformation, so we can make the “g.o.f.” statistic more or less anything at all—a serious pathology.

At this point, experts point out that *ratios* of likelihoods have the desired invariance under change of variable, but, while the likelihood ratio is a useful test statistic in certain special cases, it is not at all clear how to obtain a useful g.o.f. statistic from the likelihood ratio in the general, unbinned, case.

## V A REPLACEMENT MODEL

Since we now lack an intuitive understanding, we need a replacement intuition for what is going on. I propose this model:

Denote by $`H_0`$ the hypothesis that the data are from the p.d.f. in question. Specify an alternative hypothesis $`H_1`$ that the data are from a uniform p.d.f. (flat in the variables that we happen to have chosen). At least, the $`H_1`$ p.d.f. is flat over the region where we have data—outside that region it can be cut off.

Performing a classic Neyman-Pearson hypothesis test of $`H_0`$ vs $`H_1`$, we use the ratio of their likelihoods as our test statistic:

$$\lambda (\stackrel{}{x})=\frac{L(\stackrel{}{x}|H_0)}{L(\stackrel{}{x}|H_1)}=\frac{L(\stackrel{}{x})}{\mathrm{constant}}$$

So, the “g.o.f.” statistic can be re-interpreted as suitable for a hypothesis test that indicates which of $`H_0`$ (our p.d.f.) and $`H_1`$ (a flat p.d.f.) is more favored by the data—a well established statistical practice.

The benefit of the new interpretation is that it explains behaviors that were baffling under the g.o.f. interpretation: Neyman-Pearson hypothesis tests and g.o.f. tests behave quite differently.

For example, a reasonable g.o.f. statistic should be at least approximately distribution independent, but $`\lambda (\stackrel{}{x})`$ is often highly correlated with the m.l.e.’s (100% in our exponential case). This high correlation was confirmed in the example contributed by K. Kinoshita to the 2002 Durham Conference. Not knowing the true value of the parameters then makes it difficult, or impossible, to use $`\lambda (\stackrel{}{x})`$ as g.o.f., since we don’t know what $`\lambda (\stackrel{}{x})`$ *should* be.333Small correlations are not fatal. For example, if the P-value of g.o.f. for the observed data in a particular case ranged only between, say, 20% and 30%, for different true values within $`\pm 3\sigma `$ of the estimated value of a parameter, one would be justified in concluding “good fit” (assuming the g.o.f. statistic used had the right properties in other respects). The behavior of these correlations is natural and obvious in the hypothesis test picture: changing the parameters changes the “flatness” of the $`H_0`$ p.d.f., and $`\lambda (\stackrel{}{x})`$ reflects this.

Reference  pointed out that, with no unknown parameters, one can always transform the p.d.f. to a flat distribution. Then $`\lambda (\stackrel{}{x})`$ becomes constant independent of the data—bad news for g.o.f. In the hypothesis test picture, this becomes a comparison between two identical hypotheses, and the result is what we would expect.

## VI TEST BIAS

Take the $`H_0`$ p.d.f. to be

$$e^t(t0)$$

This distribution is fully specified—no unknown parameters. Our “g.o.f.” statistic is then

$$\mathrm{ln}L=N\widehat{t}$$

whose mean is $`\mathrm{ln}L=N`$, and variance is $`\mathrm{Var}(\mathrm{ln}L)=N`$, for an ensemble of data sets from the $`H_0`$ p.d.f. A data set with $`\widehat{t}`$ close enough to 1 will be claimed to be a good fit to the $`H_0`$ p.d.f.

But say, unknown to us, the data are really from a triangular p.d.f.:

$$1|t1|(0t2)$$

The mean and variance of $`N\widehat{t}`$ will be $`N`$ and $`N/6`$ respectively, for data from the triangular distribution. So, although the exponential and triangular p.d.f.’s are quite different, the triangular data will be more likely to pass the g.o.f. test than exponential data for which it was intended. Statisticians refer to this situation as a case of “test bias”.

We conclude that, even with no free parameters, the “g.o.f.” test is biased: there exist “impostor” p.d.f.’s that should produce bad fits, but instead pass the “g.o.f.” test with greater probability then the p.d.f. for which the test was designed. Reference  gives additional examples of this behavior.

From the hypothesis test point of view, this behavior makes sense. The exponential and triangular data have the same “distance” from the flat distribution, on the average, with the triangular data being less susceptible to fluctuations. The hypothesis test doesn’t tell us when the data are inconsistent with both $`H_0`$ and $`H_1`$.

## VII ANOTHER EXAMPLE

Here we try to find an example p.d.f. (with a free parameter) that the method in question can handle well. We use the insight provided by the hypothesis test picture. We want to keep the correlation between the free parameter and the g.o.f. statistic $`L_{\mathrm{max}}`$ to a minimum. In the hypothesis test picture, this is achieved when the “flatness” of the p.d.f. is independent of the parameter. A location parameter has this property. Additionally, we want the p.d.f. to be easily distinguishable from a flat p.d.f. So we choose the Gaussian

$$\frac{1}{\sqrt{2\pi }\sigma }e^{0.5(x\mu )^2/\sigma ^2}$$

where $`\mu `$ is unknown, but $`\sigma `$ is specified in advance. The likelihood is given by

$$\mathrm{ln}L=\underset{i=1}{\overset{N}{}}\left[\mathrm{ln}\sqrt{2\pi }+\mathrm{ln}\sigma +\frac{1}{2}\left(\frac{x_i\mu }{\sigma }\right)^2\right]$$

When $`\mu `$ and $`\sigma `$ are both unknown, their m.l.e.’s are

$$\widehat{\mu }=\frac{1}{N}\underset{i=1}{\overset{N}{}}x_i\widehat{\sigma }^2=\frac{1}{N}\underset{i=1}{\overset{N}{}}(x_i\widehat{\mu })^2$$

Using these expressions, we can rewrite the likelihood in the form $`L(\mu ,\sigma ;\widehat{\mu },\widehat{\sigma })`$:

$$\mathrm{ln}L=\frac{N}{2}\left[\mathrm{ln}(2\pi )+\mathrm{ln}(\sigma ^2)+\frac{\widehat{\sigma }^2+(\widehat{\mu }\mu )^2}{\sigma ^2}\right]$$

When only $`\mu `$ is unspecified, its m.l.e. is $`\widehat{\mu }`$ as above, and the value of the maximized likelihood is

$$\mathrm{ln}L_{\mathrm{max}}=\frac{N}{2}\left[\mathrm{ln}(2\pi )+\mathrm{ln}(\sigma ^2)+\frac{\widehat{\sigma }^2}{\sigma ^2}\right]$$

Our victory is that $`L_{\mathrm{max}}`$ only depends on $`\widehat{\sigma }`$, which is an ancillary statistic for $`\mu `$. That is, we don’t need to know the true value of $`\mu `$ in order to calculate the distribution of our g.o.f. statistic in this carefully chosen example. In fact, a convenient form for the g.o.f. statistic is

$$N\frac{\widehat{\sigma }^2}{\sigma ^2}=\underset{i=1}{\overset{N}{}}\left(\frac{x_i\widehat{\mu }}{\sigma }\right)^2$$

which is well known to have the distribution (under the null hypothesis) of a $`\chi ^2`$ with $`N1`$ degrees of freedom.

### VII.1 The Bad News

Before we declare that the method performs well in this example, there are several ugly facts to consider:

* Data that match the null hypothesis well yield $`N\widehat{\sigma }^2/\sigma ^2N`$. Much larger or much smaller values of the g.o.f. statistic imply poor g.o.f. This is in contrast to Pearson’s $`\chi ^2`$ (binned $`\chi ^2`$), for example, where smaller $`\chi ^2`$ is always better g.o.f. So we must interpret this statistic differently than how we are used to.
* The g.o.f. in this example simply reduces to a comparison between the sample variance and $`\sigma ^2`$. Any distribution with variance approximately equal to $`\sigma ^2`$ will usually generate data that “pass the test”, even distributions that look nothing like a Gaussian. This is the same kind of problem that we first saw in section III.1.
* A construction similar to that of section VI will produce “impostor” p.d.f.’s that pass the “g.o.f.” test with greater frequency than the null hypothesis. So, we have not eliminated the test bias problem.

In this example, the g.o.f. method in question will be able to flag some, but not all, of data samples that poorly match the null hypothesis. In answer to the question “Are the data from a Gaussian with unspecified mean, and variance equal to $`\sigma ^2`$?”, this g.o.f. method can only answer “No” or “Maybe”: it checks the variance part of the question, but does nothing to check the Gaussian part.

## VIII CONCLUSIONS

* This “g.o.f.” method is fatally flawed in the unbinned case. Don’t use it. Complain when you see it used.
* With fixed p.d.f.’s, the method suffers from test bias, and is not invariant with respect to change of variables. These problems persist when there are floating parameters.
* With floating parameters, the method is often circular: “g.o.f.” becomes a comparison between the measured values and the true (but unknown) values of the parameters…
* The misbehavior of this “g.o.f.” statistic is understandable when reinterpreted as the ratio between the likelihood in question and a uniform likelihood, and used to distinguish between these two specific hypotheses. Dual-hypothesis tests are not g.o.f. tests.

###### Acknowledgements.

I would like to thank Louis Lyons for several helpful discussions of the points raised here, and the organizers of the PHYSTAT2003 Conference for arranging a superb program.