text:
# Bose-Einstein condensation in a stiff TOP trap with adjustable geometry

## Abstract

We report on the realisation of a stiff magnetic trap with independently adjustable trap frequencies, $`\omega _z`$ and $`\omega _r`$, in the axial and radial directions respectively. This has been achieved by applying an axial modulation to a Time-averaged Orbiting Potential (TOP) trap. The frequency ratio of the trap, $`\omega _z/\omega _r`$, can be decreased continuously from the original TOP trap value of 2.83 down to 1.6. We have transferred a Bose-Einstein condensate (BEC) into this trap and obtained very good agreement between its observed anisotropic expansion and the hydrodynamic predictions. Our method can be extended to obtain a spherical trapping potential, which has a geometry of particular theoretical interest.

Since the first observation of BEC in a dilute alkali gas in 1995 , the spectroscopy of the excitations of these novel quantum objects has been an important tool for understanding them. A close analogy can be drawn between the discrete excited states of a condensate and the well defined energy levels of an atom. Spectroscopic measurements of atomic energy levels have guided the development of the quantum mechanical model of the atom and likewise the spectroscopy of the condensate should provide crucial information about its nature. The first experiments focused on low energy collective excitations both in a pure condensate close to absolute zero  and at finite temperature . These latter experiments led to a better understanding of the interaction between the condensate and the thermal cloud . The first observation of a transverse collective excitation was recently reported in . This also gave evidence for superfluidity in a condensate. Recent theoretical work has provided information on the behaviour of an excited condensate in specific anisotropic trap geometries. This paper presents an adjustable magnetic trap that allows us to access these geometries whilst retaining a stiff trapping potential.

In the linear hydrodynamic theory of Bose-Einstein condensation, the spectrum of low energy collective excited states depends strongly on the trap geometry (i.e. the trap stiffnesses in all three orthogonal directions). By changing the trap frequency in one specific direction, the mode energies change relative to one another and new features appear in the spectrum. For smaller condensates, the full Gross-Pitaevskii equation including the kinetic energy term is required. In this regime, the size of the condensate, as well as its shape, plays an important role in determining the excitation spectrum. If the shape of the potential is adjusted so that the second harmonic of one mode corresponds to the fundamental frequency of another, then the inherent non-linearity of the system leads to coupling between the modes. It has been predicted from the hydrodynamic limit of the Gross-Pitaevskii equation that this behaviour should be observed for a trap potentials with frequency ratios of 1.95, 1.51 and 0.68 .

Another geometry of great interest is the spherically symmetric trap, because the excitation spectrum of the condensate is simplified by degeneracy and becomes theoretically tractable . One mechanism for the decay of collective excitations in a condensate is Landau damping, which is analogous to a scattering process. A condensate excitation disappears in a collision with a thermal excitation, resulting in a thermal excitation at a higher energy. The reduced availability of states in a spherical trap, caused by degeneracy, has a major effect on the Landau damping rate. Quantitative calculations of this process are currently being carried out .

In addition to its importance for studying the collective excitations of a BEC, the spherically symmetric harmonic trap also produces unique behaviour in a thermal cloud of atoms. The breathing mode oscillation is completely undamped in this geometry , whilst in an anisotropic trap all modes are strongly damped .

There are two general types of purely magnetic trap for neutral atoms. Ioffe-Pritchard type traps  use static magnetic fields to generate a very tight radial confinement. They produce prolate condensates with aspect ratios of $`0.05`$. The aspect ratio of a condensate or thermal cloud is the ratio of the radial to axial size, and is the inverse of the trap frequency ratio, in a harmonic potential. The TOP trap , in its axially symmetric configuration, combines a quadrupole field with a rotating bias field to trap an oblate condensate, with an aspect ratio of 2.83. In this paper we demonstrate a novel method to independently relax the axial trapping frequency of the TOP trap. This new trap, with its adjustable aspect ratio, has been called a zTOP trap because of the configuration of bias fields that is used.

The original attempts to vary the aspect ratio of the TOP trap made use of the force of gravity on atoms held in a very weak magnetic potential . This method, along with its practical difficulties, is discussed briefly. In contrast, the zTOP trap combines an adjustable aspect ratio with a stiff trapping potential. The experimental setup is described and finally evidence for the successful transfer of a condensate into the zTOP trap is presented.

The standard TOP trap combines a spherical quadrupole field of radial gradient $`B_q^{}`$, with a bias field, $`𝐁_𝐭`$, rotating in the horizontal xy-plane. The expression for the trap potential usually ignores gravity. However, for heavy atoms such as 87Rb, in a weak trap, it must be included. The full potential is given by the time-average of Eq. 2, where $`r_0=B_t/B_q^{}`$.

$`U(x,y,z,t)=\mu B_q^{}`$ (2)
$`\left(|\left(x+r_0\mathrm{cos}\omega _tt\right)\widehat{𝐱}+\left(y+r_0\mathrm{sin}\omega _tt\right)\widehat{𝐲}2z\widehat{𝐳}|+{\displaystyle \frac{mgz}{\mu B_q^{}}}\right)`$
At small distances from the trap centre, $`rr_0`$, this approximates to a harmonic trapping potential. The trap oscillation frequencies and frequency ratio are

$`f_r`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}\sqrt{{\displaystyle \frac{\mu }{2m}}}{\displaystyle \frac{B_q^{}}{\sqrt{B_t}}}(1\eta ^2)^{1/4}(1+\eta ^2)^{1/2}`$ (3)
$`f_z`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}\sqrt{{\displaystyle \frac{\mu }{2m}}}{\displaystyle \frac{B_q^{}}{\sqrt{B_t}}}\sqrt{8}(1\eta ^2)^{3/4}`$ (4)
$`{\displaystyle \frac{f_z}{f_r}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{8(1\eta ^2)}{1+\eta ^2}}}`$ (5)
where the parameter

$$\eta =\frac{mg}{2\mu B_q^{}}$$
(6)
is the ratio of the gravitational force on the atom to the axial magnetic force due to the quadrupole field. Note that the trap frequency ratio depends only on $`\eta `$. Thus by measuring the axial and radial trap frequencies, as a function of $`B_q^{}`$, we were able to accurately calibrate the quadrupole field strength.

However, the aspect ratio only changes significantly in very weak traps, where the gravitational force becomes important. A ten percent decrease in the aspect ratio requires $`B_q^{}<23`$ G/cm. Since $`B_t`$ is limited by noise to a minimum value of 2 G in our apparatus, this corresponds to a radial trapping frequency of $`<15`$ Hz. Attempts were made to selectively excite quadrupole oscillations of the BEC in this weak gravity trap. However, in such a weak trap it was impossible not to excite large dipole motions, making the quadrupole oscillation impossible to observe.

The theory of the zTOP trap is analogous to that for the standard TOP trap, except that the bias field oscillates in three rather than two dimensions. Both are based on a spherical quadrupole field trap, with a trapping potential of

$$U(x,y,z)=\mu B_q^{}\left|x\widehat{𝐱}+y\widehat{𝐲}2z\widehat{𝐳}\right|$$
(7)
which has a restoring force in the axial direction that is twice as strong as that in the radial direction. Now consider the effect of applying a TOP bias field rotating at $`\omega _t`$ in the (radial) xy-plane. This causes the locus of the quadrupole field to describe a circle in the xy-plane and after time-averaging creates an axially symmetric TOP trap. The trap stiffness is reduced in all directions compared to the quadrupole field trap, but more significantly in the radial direction, increasing the aspect ratio to 2.83. Extending this argument, the trap stiffness can be preferentially reduced in the z direction by including an oscillating axial bias field. Our present arrangement uses an axial modulation frequency ($`\omega _a`$) of $`2\omega _t`$. Thus in our zTOP trap, the locus of the quadrupole field follows a three dimensional saddle shape rather than a circle. The total zTOP magnetic field is given in Eq. 9 as the sum of its three components - the static quadrupole field and the oscillating radial and axial bias fields.

$`𝐁(t)=B_q^{}\left(x\widehat{𝐱}+y\widehat{𝐲}2z\widehat{𝐳}\right)+`$ (9)
$`B_t\left(\mathrm{cos}\omega _tt\widehat{𝐱}+\mathrm{sin}\omega _tt\widehat{𝐲}\right)+B_z\mathrm{cos}\omega _at\widehat{𝐳}`$
The properties of the zTOP trap are calculated from the instantaneous zTOP potential, numerically averaged over one cycle of the bias field. This calculation gives the following constraints on the frequency, $`\omega _a`$, of the axial modulation necessary to create a stable trap. Firstly, to create a time-averaged potential, $`\omega _a`$ must be greater than the trap frequencies. Secondly $`\omega _a`$ must be an exact integer multiple, n, of the radial (TOP) bias field frequency, $`\omega _t`$. Our calculations showed that deviations from this condition produced asymmetry in the time-averaged potential between the x and y directions. In the xy-plane, the trap will acquire a slightly elliptical, rotating cross-section, whilst its centre will describe a small circle. These slow variations in the time-averaged potential at the beat frequency, $`\omega _a`$ \- n$`\omega _t`$, could heat the cloud. It is also possible to choose $`\omega _t`$ as an integer multiple of $`\omega _a`$. However this configuration requires the two signals to be phase locked, if the ’micromotion’ of the trap described above is to be avoided. In addition, $`\omega _a`$ must be lower than the Larmor frequency, so that Majorana spin flips are not induced.

The extent of the trap relaxation in the axial direction depends on the waveform and amplitude, $`B_z`$, of the axial modulation, but not on its frequency. The axial modulation causes the locus of the quadrupole field to oscillate in the vertical direction. A waveform that causes it to spend more time per cycle at the extremes of this motion results in a more relaxed trap in the axial direction. Thus an ideal square wave produces the weakest trap, for a given maximum $`B_z`$ field. Removing the higher harmonics in the axial modulation stiffens the trap as shown in Fig.1. Note that significant additional relaxation, with respect to single frequency modulation, is achieved by adding only the third harmonic component of a square wave to the fundamental. In this case a spherical trap is produced for $`B_z/B_t=5.5`$. Since $`B_t`$ is typically 2 G, this is a readily achievable experimental condition. The axial trapping frequency is also reduced by increasing the modulation amplitude, $`B_z`$. Thus we are able to vary the aspect ratio of the trap during the course of an experiment by controlling the peak voltage across the axial bias coils.

The apparatus that we use to create 87Rb condensates is described in detail elsewhere . In summary, we use a differentially pumped double MOT system (which incorporates a pyramidal MOT ) to load a TOP trap with $`2\times 10^8`$ atoms. Evaporative cooling proceeds via both Majorana spin flip and radio frequency cutting to achieve condensates with $`2\times 10^4`$ atoms.

To generate the axial bias field, two Helmholtz coils of 35 turns each, have been added above and below the experimental cell. A 2 turn pickup coil is also present to monitor the applied field. We currently use an axial modulation at 14 kHz (2$`\omega _t`$). The amplitude of this signal is computer controlled and fed to an audio amplifier. A transformer is used to match the output of the amplifier to the coils. Fine adjustments are made to the frequency until the x and z pickup coils produce a clean Lissajous figure on the oscilloscope, indicating that the axial oscillation is at exactly twice the frequency of the radial bias field. This frequency relationship is stable to within $`1`$ Hz over the day.

The trap frequencies were calibrated with respect to the voltage supplied across the axial bias field coils. A small, cold thermal cloud was collected in a ’displaced’ zTOP trap. The additional static magnetic fields causing the displacement were suddenly switched off, exciting orthogonal dipole oscillations at the two trap frequencies. The radial and axial positions of the cloud were recorded as a function of time, using absorption imaging, and fitted with sine waves to obtain the trap frequencies. Figure 2 shows the data for the normalised trap frequencies and the trap frequency ratio, for a range of axial bias field amplitudes.

Theoretical values for the normalised trap frequencies, as a function of $`B_z/B_t`$, were found by numerically integrating the magnetic potential energy of an atom over one cycle of the bias field. These values were accurately described by an eighth order polynomial. Two free fitting parameters were used to fit this polynomial to the experimental data. One parameter determined the trap frequencies without any axial bias field and the other was a constant, relating the voltage across the axial bias coils to the actual value of $`B_z/B_t`$ that was produced. Fitting this expression to the axial trap frequency data gave values for both fitting parameters and their respective errors. The solid lines in Fig.2 show the fitted zTOP trap calibration curves.

The condensate is initially produced by a sequence of computer controlled evaporative cooling ramps in a standard TOP trap, with final trap frequencies $`\omega _r=126`$ Hz and $`\omega _z=356`$ Hz $`(=\sqrt{8}\omega _r)`$. The amplitude of the oscillating axial field is then linearly ramped from zero to its final value over 0.5s. This adiabatically transfers the condensate into the zTOP trap.

Figure 3 shows the aspect ratio of a BEC, 12ms after release from zTOP traps with a range of frequency ratios. Changing the frequency ratio of the trap changes the distribution of the self energy of the condensate among the three orthogonal directions and hence the form of the free expansion of the condensate. The aspect ratio of the expanding BEC has been predicted as a function of trap geometry using the hydrodynamic equations  and is shown as a solid line in Fig. 3. This theoretical curve contains no free parameters. The good agreement between theory and experiment confirms two facts. Firstly that the condensed state survived adiabatic transfer to the zTOP trap and secondly that the trap potential experienced by the condensate has been modified as predicted.

In conclusion, we have modified our TOP trap to create a stiff zTOP trap with an adjustable geometry. Trap frequency ratios between 2.83 and 1.6 have been demonstrated by applying a magnetic field in the axial direction, oscillating at a single frequency. The addition of higher harmonics to the axial bias field modulation will make a spherical trap accessible. This prediction, combined with the successful transfer of BECs into the zTOP trap, as demonstrated in this paper, promises to open up the field of degenerate condensate excitations to experimental reasearch.

This work was supported by the EPSRC and the TMR program (No. ERB FMRX-CT96-0002). O.M. Maragò acknowledges the support of a Marie Curie Fellowship, TMR program (No. ERB FMBI-CT98-3077).


format: text

text:
# An extreme critical space-time: echoing and black-hole perturbations

## I Introduction

The discovery by Choptuik of critical phenomena in gravitational collapse is widely regarded as the most unexpected success of numerical relativity. The now famous mass scaling law and echoing behavior had not been anticipated by existing black-hole theory. Despite a subsequent theoretical framework involving dynamical-systems and renormalization-group methods, this is still not fully understood; see reviews by Choptuik and Gundlach. In particular, the exactly critical space-time suggested by such numerical simulations is not known analytically. Thus it is of some interest to study solutions with similar properties.

The numerical results suggest a self-similar critical space-time. For the massless Einstein-Klein-Gordon case originally studied by Choptuik, there is a discrete self-similarity with a characteristic echoing period. Other matter models yield continuous self-similarity, which can be described by a homothetic vector. A class of homothetic, spherically symmetric solutions to the massless Einstein-Klein-Gordon equations was given by Roberts and has been studied by various authors. This article describes another such solution, which is static. This solution, the Roberts class and its time-inverse form a two-parameter family which was actually given by Brady, though in that reference a parameter was fixed to reduce to the Roberts class.

Since these solutions are continuously rather than discretely self-similar, they do not include the critical space-time itself. However, one may hope to find a perturbed space-time close to the critical space-time. Indeed, Frolov has argued that linear perturbations of the critical Roberts solution break the self-similarity from continuous to discrete. This article therefore studies the static solution and its spherically symmetric linear perturbations. There are indeed perturbations which break the continuous self-similarity to a certain periodic structure, actually falling short of discrete self-similarity for reasons explained subsequently. The static symmetry greatly simplifies the analysis. In particular, this allows the first explicit construction of perturbed black holes with finite mass from such a critically trapped space-time.

The article is organised as follows. Section II gives the solution and describes the global structure of the space-time. Section III studies linear perturbations, explicitly giving the general solution to the linearized field equations, then imposing appropriate boundary conditions and finding the relevant modes. Section IV describes a type of dominant mode which produces echoing. Section V describes the other dominant mode, which produces black holes for a certain quadrant in amplitude space. Section VI shows that the other quadrants produce white holes, naked singularities and regular dispersal respectively. Section VII concludes with some speculations concerning critical collapse.

## II The solution

The field equations for Einstein gravity with a massless Klein-Gordon field $`\varphi `$, fixing the relative coupling constant, are

$`R=2\varphi \varphi `$ (1)
$`^2\varphi =0`$ (2)
where $`R`$ is the Ricci tensor and $``$ the covariant derivative of the metric $`g`$. A spherically symmetric metric may be written locally in terms of the line-element

$$ds^2=r^2d\mathrm{\Omega }^22e^{2\gamma }dx^+dx^{}$$
(3)
where $`d\mathrm{\Omega }^2`$ is the line-element of the unit sphere and $`(r,\gamma )`$ are functions of the null coordinates $`x^\pm `$. Then the area of the spheres is $`4\pi r^2`$, so that $`r`$ is the areal radius. The field equations in these coordinates are

$`_\pm _\pm r2_\pm \gamma _\pm r=r(_\pm \varphi )^2`$ (4)
$`r_+_{}r+_+r_{}r+e^{2\gamma }/2=0`$ (5)
$`r^2_+_{}\gamma _+r_{}re^{2\gamma }/2=r^2_+\varphi _{}\varphi `$ (6)
$`r_+_{}\varphi +_+r_{}\varphi +_{}r_+\varphi =0.`$ (7)
It will also be convenient to use coordinates $`(\rho ,\tau )`$ defined locally by

$$x^\pm =\pm e^{\rho \pm \tau }.$$
(8)
The solution is given simply by

$$r=e^\rho \gamma =0\varphi =\tau $$
(9)
which is the case $`\alpha =\beta =0`$ of Brady. The variables $`(X,Y)`$ of Choptuik evaluate as $`(0,1/2)`$. The line-element can be written as

$$ds^2=e^{2\rho }(d\mathrm{\Omega }^2+2d\rho ^22d\tau ^2).$$
(10)
This is explicitly conformal to a flat metric times a unit sphere, so that a Penrose diagram can be given as in Fig.1. The metric is evidently static, with static Killing vector $`/\tau `$. It is also homothetic, or continuously self-similar, with homothetic vector $`/\rho `$. This means that

$$L_{/\rho }g=2g$$
(11)
where $`L`$ denotes the Lie derivative and the factor of 2 is conventional. It follows that the null vectors $`/\xi ^\pm `$ are also homothetic vectors, where

$$\xi ^\pm =\pm \mathrm{ln}(\pm x^\pm )=\tau \pm \rho .$$
(12)
In this sense the solution is doubly self-similar. Hypersurfaces orthogonal to a homothetic vector may be called homothetic hypersurfaces.

There is a central singularity, meaning $`r=0`$ there, at $`x^+x^{}=0`$. This is a bifurcate null singularity, where the causal nature is defined by the conformal metric. It is therefore also marginally trapped, since $`_\pm r=0`$ on $`x^{}=0`$. It is a curvature singularity, since the Ricci scalar is found to be just $`r^2`$. The other conformal boundary corresponds to infinite $`r`$ and constitutes null infinity, though the space-time is not asymptotically flat. Correspondingly, the mass or energy

$$m=(r/2)(1g^1(r,r))$$
(13)
evaluates as $`m=r/4`$, becoming infinite at null infinity. The singularity is massless in the sense that $`m=0`$ there. This confirms the view that massless central singularities lie at the threshold between generic trapped and naked singularities, which occur for $`m>0`$ and $`m<0`$ respectively. The space-time is summarized in Fig.1. The $`(\rho ,\tau )`$ or $`\xi ^\pm `$ coordinates take all real values, whereas the $`x^\pm `$ coordinates have the range $`x^+>0`$, $`x^{}<0`$. The future and past branches of the singularity occur for $`x^{}=0`$ and $`x^+=0`$ respectively, where $`/\tau `$ is taken to be future-pointing.

The solution just avoids being a black hole (or white hole) in the sense that there are no trapped surfaces, yet the singularity is marginally trapped. Also, the singularity just avoids being naked in the sense of Penrose, that it should lie to the future of some point and the past of another point on some observer’s worldline. Likewise, an observer could see arbitrarily high but not infinite curvature from the future branch of the singularity. In summary, this is a new type of extreme space-time, lying at the threshold between black holes and naked singularities.

## III Linear perturbations

Henceforth, spherically symmetric linear perturbations of the space-time will be studied, initially following the analysis of the critical Roberts solution by Frolov. The linear perturbations studied here may be decomposed by Laplace transformations into modes

$`r=e^\rho (1+ϵ\stackrel{~}{r}(\tau )e^{k\rho })`$ (14)
$`\gamma =ϵ\stackrel{~}{\gamma }(\tau )e^{k\rho }`$ (15)
$`\varphi =\tau +ϵ\stackrel{~}{\varphi }(\tau )e^{k\rho }`$ (16)
governed by the complex frequency $`k`$. Assuming this form of the perturbations fixes the perturbative gauge freedom, apart from the zero point of $`\tau `$. The field equations are linearized by removing terms $`o(ϵ)`$ in the perturbation parameter, $`0<ϵ1`$. The four independent components of the field equations yield three independent linearized equations:

$`\stackrel{~}{\gamma }=(1k/2)^2\stackrel{~}{r}\stackrel{~}{r}^{\prime \prime }/4`$ (17)
$`\stackrel{~}{\varphi }=\stackrel{~}{r}^{}\stackrel{~}{\gamma }^{}/k`$ (18)
$`\stackrel{~}{r}^{\prime \prime }+k^2\stackrel{~}{r}+2k\stackrel{~}{\gamma }+2\stackrel{~}{\varphi }^{}=0`$ (19)
where $`f^{}=f/\tau `$. The first two are solved for $`(\stackrel{~}{\gamma },\stackrel{~}{\varphi })`$ explicitly if $`\stackrel{~}{r}`$ is known, so that the system reduces to a single fourth-order equation for $`\stackrel{~}{r}`$:

$$\stackrel{~}{r}^{\prime \prime \prime \prime }2(k^2k+2)\stackrel{~}{r}^{\prime \prime }+k^2(k^22k+4)=0.$$
(20)
This linear ordinary differential equation can be solved via the auxiliary equation

$$\omega ^42(k^2k+2)\omega ^2+k^2(k^22k+4)=0.$$
(21)
This is a biquadratic equation, so that the roots come in pairs, $`\omega =\pm \omega _A`$, $`\pm \omega _B`$ where

$`\omega _A=(k^22k+4)^{1/2}`$ (22)
$`\omega _B=k.`$ (23)
Then the general solution for each $`k`$ is

$$\stackrel{~}{r}=A_+e^{\omega _A\tau }+A_{}e^{\omega _A\tau }+B_+e^{\omega _B\tau }+B_{}e^{\omega _B\tau }$$
(24)
for constants $`A_\pm `$, $`B_\pm `$, except in special cases where roots coincide, when similarly standard solutions can be given. The special cases are (i) $`k=0`$, when $`\omega _B=0`$, (ii) $`k=2`$, when $`\omega _A=\omega _B`$, and (iii) $`k=1\pm \sqrt{3}i`$, when $`\omega _A=0`$.

Weak boundary conditions will be taken as

$$\underset{x^\pm \pm \mathrm{}}{lim}r0\underset{x^\pm 0}{lim}r^10.$$
(25)
If these conditions were not met, the perturbed $`r`$ would vanish at the original infinity or become infinite at the original centre, respectively, either of which would be wildly different from the background space-time and therefore inconsistent with the linear approximation. Slightly stronger boundary conditions additionally impose

$$\underset{x^\pm \pm \mathrm{}}{lim}r^1\underset{x^\pm 0}{lim}r$$
(26)
with similar justification. By writing a mode explicitly as

$`r`$ $`=`$ $`e^\rho (1+ϵAe^{\omega \tau k\rho })`$ (27)
$`=`$ $`(x^+x^{})^{1/2}+ϵA(x^+)^{(1k+\omega )/2}(x^{})^{(1k\omega )/2}`$ (28)
the weak boundary conditions can be seen to imply

$$\mathrm{}(k\pm \omega )1$$
(29)
with cases of equality examined individually in the following. For the $`\omega _B`$ modes this means

$$a1/2$$
(30)
where $`a=\mathrm{}k`$, $`b=\mathrm{}k`$, whereas for the $`\omega _A`$ modes a little calculation yields

$$a=1b^23$$
(31)
with $`\omega _A=i\beta `$, $`\beta ^2=b^23`$. These modes are indicated in Fig.2.

## IV Echoing

A common argument is that the relevant modes are those with $`a>0`$ and the dominant modes those with maximal $`a`$. As depicted in Fig.2, this occurs for the $`\omega _A`$ modes. One might also expect the dominant mode among these to be that with minimal $`|b|`$, since higher-frequency modes interfere destructively. This mode is given by $`k=1\pm \sqrt{3}i`$, case (iii) above:

$$r=e^\rho +ϵA\mathrm{cos}(\sqrt{3}\rho +c).$$
(32)
The general solution for this mode would replace $`A`$ with $`A+C\tau `$, but if $`C0`$ the boundary conditions (25) are not satisfied. The perturbed space-time is static and the perturbation $`\delta r`$ is periodic, with echoing in $`\rho `$ of period

$$\mathrm{\Delta }=2\pi /\sqrt{3}=3.62\mathrm{}$$
(33)
which is within about $`5\%`$ of the value 3.44 obtained numerically by Choptuik. This suggests that the perturbed space-time may be close to the critical space-time. However, the Klein-Gordon perturbation $`\delta \varphi `$ vanishes for this mode.

Note that the perturbed space-time is not discretely self-similar; despite $`\delta r`$ being periodic in $`\rho `$, $`r`$ contains an exponential term in $`\rho `$. The situation may be understood schematically in terms of a phase-space picture. The universality of the critical exponent suggests that the critical solution is an intermediate attractor of codimension one. In other words, there is a critical hypersurface in phase space, separating black-hole solutions from dispersal solutions. The background solution lies in this hypersurface, by the same argument as for the critical Roberts solution: the Roberts class, cut and pasted to flat regions as described by Brady, describe either black-hole formation or dispersal, except for the delimiting critical case. Linear perturbations are represented by vectors from the background solution. One seeks a perturbation tangent to the critical hypersurface, pointing towards the critical space-time. If the critical hypersurface is curved, the perturbed solution will generally be outside the hypersurface and so exhibit exponential as well as oscillatory behavior.

For the modes with $`k=1+ib`$, $`b^2>3`$, one finds

$$r=e^\rho +ϵA\mathrm{cos}(b\rho +c)\mathrm{cos}(\beta \tau +d)$$
(34)
so that the perturbations are periodic in both $`\rho `$ and $`\tau `$. The variable $`X`$ plotted by Choptuik is, to $`O(ϵ)`$,

$$\frac{1}{2}\frac{\varphi }{\rho }=\frac{ϵAb\beta }{4}\mathrm{sin}(b\rho +c)\mathrm{sin}(\beta \tau +d)$$
(35)
which is explicitly sinusoidal in both $`\rho `$ and $`\tau `$. Apart from the well known echoing in $`\tau \rho `$, roughly sinusoidal oscillations in $`\rho `$ at constant $`\tau `$ do seem to develop in the numerical simulations, e.g. figure 3 of Choptuik.

## V Black holes

The common argument that the dominant mode will occur for maximal $`\mathrm{}k`$, as this is largest as $`\rho \mathrm{}`$, implicitly assumes that the other coordinate $`\tau `$ is finite, which actually applies only where the singularity bifurcates. For the rest of the singularity, the explicit expression (28) shows that the dominant mode at $`x^{}=0`$ instead occurs for maximal $`\mathrm{}(k\pm \omega )`$. As above, this maximal value is unity. The corresponding modes consist of those discussed in the previous section and the $`\omega _B`$ modes with $`\mathrm{}k=1/2`$. Non-linear perturbations may then be dominated by either or both, in a way which is unclear from a linear perturbation analysis. Henceforth the new modes will be discussed separately.

The $`\omega _B`$ modes with $`k=1/2+ib`$, $`b0`$, do not satisfy the boundary conditions (26), since $`r`$ does not have a limit as $`x^\pm 0`$, $`x^{}0`$. Thus the dominant mode is just $`\omega _B=k=1/2`$. Then

$`r`$ $`=`$ $`e^\rho +ϵB_+e^{(\rho +\tau )/2}+ϵB_{}e^{(\rho \tau )/2}`$ (36)
$`=`$ $`(x^+x^{})^{1/2}+ϵB_+(x^+)^{1/2}+ϵB_{}(x^{})^{1/2}.`$ (37)
The locus of the centre $`r=0`$ can be determined by noting that it is a hyperbola in $`z^\pm =(\pm x^\pm )^{1/2}`$ coordinates:

$$0=z^+z^{}+ϵB_+z^++ϵB_{}z^{}.$$
(38)
Likewise, the locus of the trapping horizon, if any, is determined by

$$_\pm r=\pm (\pm x^\pm )^{1/2}((x^{})^{1/2}+ϵB_\pm )/2$$
(39)
where $`_\pm =/x^\pm `$. Here a trapping horizon is a hypersurface foliated by marginal surfaces, $`_+r=0`$ or $`_{}r=0`$, which respectively characterize black or white holes. These loci can be seen to be straight lines of constant $`x^{}`$ or $`x^+`$ respectively, or non-existent, depending on the signs of $`B_\pm `$. These lines correspond to the asymptotes

$$z^{}=ϵB_\pm $$
(40)
of the above hyperbola, which passes through $`z^+=z^{}=0`$. The part of the hyperbola and its asymptotes inside $`z^+>0`$, $`z^{}>0`$ can be mapped back to the $`x^\pm `$ coordinates, the loci depending qualitatively only on the signs of $`B_\pm `$. This, combined with the nature of the background conformal boundaries, determines the qualitative global structure of the perturbed space-time.

For the quadrant $`B_+<0`$, $`B_{}>0`$, assumed in the remainder of this section, there is a black hole, with a region of trapped surfaces following a trapping horizon and preceding a spatial central singularity, as depicted in Fig.3. A mode with $`b=0`$, $`0<a<1/2`$ is also depicted for comparison, the analysis being similar. In that case, the trapping horizon is spatial, appears at zero area and evolves towards infinite area at conformal infinity. However, for the dominant mode, a remarkable feature is that the trapping horizon is null, occuring at constant $`x^{}`$, with constant areal radius

$$r=ϵ^2B_+B_{}.$$
(41)
This seems to be an example of an isolated horizon which is not a Killing horizon; the perturbed space-time is not static. It is also clearly an event horizon, generalized from the usual definition for asymptotically flat space-times.

A curious feature is that the past branch of the background singularity is regularized by this perturbation: the $`x^\pm `$ coordinates break down at $`x^+=0`$, but it can be shown to have finite area, finite mass $`m=r/2`$ and vanishing Ricci scalar. (Essentially, $`\gamma \mathrm{}`$ there, so norms $`ff=2e^{2\gamma }_+f_{}f`$ tend to vanish, the exponential term suppressing negative powers). This implies that it is a degenerate trapping horizon. Extension inside this past event horizon seems, however, beyond the validity of the linear approximation; it is meaningless to have a linear perturbation of a region which does not exist in the background space-time.

Since the mass or energy (13) is $`m=r/2`$ on a trapping horizon, there is a mass scaling law

$$m=ϵ^2B_+B_{}/2.$$
(42)
Thus the mass is explicitly linear in the perturbation amplitudes $`|ϵB_\pm |`$ of both ingoing and outgoing modes, so that the critical exponent is unity. The renormalization-group argument, that the critical exponent is the reciprocal of the maximal $`\mathrm{}k`$, would not give this value unless $`\mathrm{}k`$ is replaced by $`\mathrm{}(k+\omega )`$, in accordance with the argument at the beginning of this section. This reflects the double self-similarity of the solution: the relevant homothetic vector at $`x^{}=0`$ is $`\pm /\xi ^\pm `$, rather than $`/\rho `$.

## VI Dispersal, naked singularities and white holes

The other quadrants of the dominant $`\omega _B`$ mode, $`k=1/2`$, may be similarly analysed. The quadrant $`B_+>0`$, $`B_{}<0`$ is a time-reverse of the black-hole case, yielding a white hole. For the quadrant $`B_+<0`$, $`B_{}<0`$, there is an eternally naked central singularity, with no trapped region outside the singularity. For the quadrant $`B_+>0`$, $`B_{}>0`$, both branches of the background singularity are regularized: the area vanishes only at $`x^+=x^{}=0`$, $`x^\pm =0`$ otherwise having finite area, finite mass $`m=r/2`$ and vanishing Ricci scalar, by a similar argument as for the black-hole case. Extension beyond these horizons again seems beyond the validity of the linear approximation. As far as can be judged from the unextended region, there is no singularity or trapped region, with the perturbed solution describing regular dispersal. The dispersal and naked-singularity solutions are depicted in Fig.4, the corresponding white-hole diagram being a time-reverse of Fig.3(ii).

This indicates that, in the phase-space picture, the background space-time lies not only in the critical hypersurface separating black-hole formation from dispersal, but also in another hypersurface separating black-hole formation from naked singularities, indicating a fourth region presumably describing white-hole decay. The picture would be analogous to Fig.5, with the critical space-time presumably lying roughly in the direction analogous to a third axis $`A`$ for the dominant $`\omega _A`$ mode.

## VII Remarks

The solution studied here describes a new type of extreme space-time, just on the verge of containing a black hole, white hole or naked singularity, without actually containing any of them. Each can be explicitly constructed by linear perturbations, as can regular dispersal solutions, so that the space-time is also critical; it lies at the threshold between black holes, white holes, naked singularities and dispersal. The symmetry and simplicity of the solution suggest it as a standard example of an extreme critical space-time. As such, it may be useful in the context of cosmic censorship. It should be noted that none of the perturbations studied here indicates a violation of cosmic censorship; the naked-singularity perturbations involve an eternally naked singularity, not one forming from regular initial conditions.

Various properties of the solution or its perturbations are analogous to those of near-critical numerical solutions in gravitational collapse. In particular, the dominant black-hole perturbation produces a black hole of constant mass and area, explicitly related to the perturbation amplitudes by a scaling law. Other perturbations produce echoing reminiscent of the numerical solutions.

It may be conjectured that the numerically suggested critical space-time has a similar global structure to that of the solution described here, though this contradicts the analysis of Gundlach. Of course, the global structure of a near-critical numerical solution is quite different, asymptotically flat with an initially regular centre. The critical and near-critical space-times should agree only close to the future branch of the singularity. A null central singularity in spherical symmetry is marginally trapped and therefore lies at the threshold of black-hole or white-hole formation, for the future-null and past-null cases respectively. An observer close to a future-null singularity could see arbitrarily high but not infinite curvature, consistently with the numerical results. Likewise, a future-null singularity is not naked, in the sense that there is no future-causal curve from it, but is just on the verge of nakedness; to coin a phrase, a teasing singularity.

Moreover, Christodoulou showed analytically that, for the massless Einstein-Klein-Gordon system, there are non-generic collapses to a central null singularity and non-generic black holes preceded by a central null singularity, as depicted in Fig.6. This suggests the following picture of critical collapse: as a family of black-hole solutions approaches criticality, the trapped region is zipped up between the trapping horizon and the centre. The closed zip is both a centre $`r=0`$ and, in a conformal sense, a trapping horizon $`_+r=0`$, and therefore is a central null singularity.

Acknowledgements. Thanks to Matt Choptuik, Andrei Frolov and Carsten Gundlach for discussions and Abhay Ashtekar and the Center for Gravitational Physics and Geometry for hospitality. Research supported by the National Science Foundation under award PHY-9800973.


text:
# References

Single particle spectrum and binding energy of nuclear matter

M. Baldo$``$ and A. Fiasconaro$``$

$``$ INFN, Sez. Catania, Corso Italia 57 - 95129 Catania, Italy

$``$ Dipartimento di Fisica, Universita’ di Catania, Corso Italia 57, 95129 Catania, Italy

Abstract. In non-relativistic Brueckner calculations of nuclear matter, the self-consistent single particle potential is strongly momentum dependent. To simplify the calculations, a parabolic approximation is often used in the literature. The variation in the binding energy value introduced by the parabolic approximation is quantitatively analyzed in detail. It is found that the approximation can introduce an uncertainty of 1-2 MeV near the saturation density.

PACS: 21.65.+f , 24.10.Cn , 45.50.Jf , 21.30.-x , 26.60.+c

1. Introduction It is one of the fundamental issue in nuclear physics to evaluate the nuclear matter binding energy and saturation properties, starting from a realistic nucleon-nucleon (NN) interaction with no free parameter. This old project requires the solution of a complex many-body problem, and has received several contributions and improvements along the years, beginning as far back as the middle of the last century. One of the main approaches to this long standing problem is the so-called hole-expansion or Bethe-Brueckner-Goldstone (BBG) theory . The first real breakthrough in this scheme was the introduction of the self-consistent single particle potential at the two hole-line level of approximation, which is then usually referred as the Brueckner-Hartree-Fock (BHF) approximation . The introduction of the self-consistent potential drastically improves the results. In particular, the binding energy and saturation density, which otherwise would turn out unreasonable, move to values which can be considered an acceptable starting approximation. The remaining discrepancy could be summarized in the celebrated “Coester band” , along which the saturation points for different NN interactions were clustering and which misses the empirical region ( corresponding approximately to a binding energy per nucleon of -16 MeV and a nucleon density of 0.17 $`fm^3`$ ). Later, the Liége group stressed  the relevance of the choice of the single particle potential. In particular they suggested the use of the “continuous choice”, which indeed appears to move the saturation point towards the empirical one, but still missing it . A period of major developments took place in the latest two decades. Starting from the works by B.D. Day , the hole-line expansion was analysed up to the three hole-line level of approximation. A strong indication of convergence of the expansion was obtained . Furthermore, BHF calculations with the continuous choice seem to get a substantially smaller corrections from three-body correlations . The results confirm that the empirical saturation point is still missed, and therefore that three-body forces are needed in the nuclear hamiltonian . In the meanwhile the relativistic Dirac-Brueckner (DB) method was developed , which already at the two hole-line level of approximation appears to be able to reproduce the empirical saturation point. The main relativistic correction introduced by the DB method is due to the structure of the Dirac 4-spinors, which in the medium appear “rotated” with respect to the free ones. The non-relativistic three-body forces and this relativistic effect of the DB approach are probably two faces of the same dynamical effect . The many-body theory has reached, therefore, such a precision that it is possible to test the nuclear hamiltonian. Because of that, the time seems to be appropriate to check the reliability of the approximations which are commonly employed in BBG calculations of nuclear matter.

In this letter we consider the BHF in the continuous choice and we analyze quantitatively the uncertainty of the results which comes out by approximating the single particle self-consistent potential with a parabolic form. This approximation is quite popular, since it allows to calculate the potential, at each iteration, only for few momenta, thus reducing drastically the computer time. The single particle potential, as obtained from fully self-consistent BHF calculations, is indeed strongly momentum dependent and not necessarily so simple as a parabola.

2. Sketch of the formalism In the BHF approximation, the nuclear matter total energy $`E`$ is obtained from the Brueckner G-matrix $`G(\omega )`$ according to the equation

$$E=\underset{k_1<k_F}{}\frac{\mathrm{}^2k_1^2}{2m}+\frac{1}{2}\underset{k_1,k_2<k_F}{}k_1k_2|G(e_{k_1}+e_{k_2})|k_3k_4_A$$
(1)
with $`|k_1k_2_A=|k_1k_2|k_2k_1`$. Here $`k_F`$ is the Fermi momentum, the summation over the momenta $`k_i`$ include spin and isospin variables. The single particle energies $`e_k`$, appearing in the entry energy of the G-matrix, are given by

$$e(k)=\frac{\mathrm{}^2k^2}{2m}+U(k)$$
(2)
where the single particle potential $`U(k)`$ is determined by the self-consistent equation

$$U(k)=\underset{k^{}<k_F}{}kk^{}|G(e_{k_1}+e_{k_2})|kk^{}$$
(3)
The self-consistency is coupled with the integral equation for the G-matrix

$$\begin{array}{cc}\hfill k_1k_2|G(\omega )|k_3k_4& =k_1k_2|v|k_3k_4+\hfill \\ & \\ \hfill +_{k_3^{}k_4^{}}k_1k_2|v|k_3^{}k_4^{}& \frac{\left(1\mathrm{\Theta }_F(k_3^{})\right)\left(1\mathrm{\Theta }_F(k_4^{})\right)}{\omega e_{k_3^{}}e_{k_4^{}}}k_3^{}k_4^{}|G(\omega )|k_3k_4\hfill \end{array}$$
(4)
where $`\mathrm{\Theta }_F(k)=1`$ for $`k<k_F`$ and is zero otherwise. The product $`Q(k,k^{})=(1\mathrm{\Theta }_F(k))(1\mathrm{\Theta }_F(k^{}))`$, appearing in the kernel of Eq. (4), enforces the scattered momenta to lie outside the Fermi sphere, and it is commonly referred as the “Pauli operator”. The self-consistent set of equations are usually solved by an iteration procedure. The G-matrix can be expanded in partial waves, according to the classification of two-nucleon channels . To avoid coupling between different two-body channels, the Pauli operator $`Q`$, as well as the two-body energies $`e_{k_3^{}}+e_{k_4^{}}`$ in the denominator, are averaged over the angle between the relative momentum $`q=(k_3^{}k_4^{})/2`$ and the total momentum $`P=k_3^{}+k_4^{}`$. Despite this approximation, which has been tested recently in ref. , the numerical solution of the coupled equations (3),(4) is quite time consuming, since the single particle potential $`U(k)`$ must be calculated in a wide range of momenta with a fine enough grid. If one assumes that the potential $`U(k)`$, or equivalently the single particle energy $`e(k)`$, has approximately a quadratic form

$$e(k)e_0+\frac{\mathrm{}^2k^2}{2m^{}}$$
(5)
then one can calculate the potential, at each iteration step, in few points only and interpolate the obtained values with a parabola. The approximation of Eq. (5) is usually called the effective mass approximation, since then the spectrum has the same shape as the free one but with an “effective mass” $`m^{}`$.

In order to test this approximation, we have performed a set of BHF calculations fully self-consistently without any assumption about the potential shape, as well as by forcing the potential to a parabolic shape by means a fitting procedure, and then compared the results.

3. Results and discussion The performed BHF calculations include all two-body channels up to total angular momentum $`J=11\mathrm{}`$. In a set of calculations we adopted the Argonne v18  potential as the NN interaction. This potential belong to a new generation of realistic NN potentials, with an improved fit of the scattering data, which give similar results and cluster closely together in the Coester band . The self-consistent single particle potential $`U(k)`$ was calculated up to the momentum cut-off $`k_{max}=7.5fm^1`$, which turns out to be large enough in the considered density range . The potential, for the Fermi momentum $`k_F=1.4fm^1`$, is displayed in Fig. 1 (full circles). It is numerically calculated with a grid step of $`0.1fm^1`$ from the G-matrix, Eq. (3), and inserted as the entry potential at each iteration step, until convergence is reached, i.e. the potential and the binding energy are stable under iteration with good accuracy. Stability within few KeV of the binding was systematically reached. The numerical method is described in ref.

The quadratic approximation, at each step of the iteration procedure, is introduced by fitting the potential up to a certain maximum momentum $`k_{FIT}`$. For definiteness, we have considered in detail two choices, namely $`k_{FIT}=2k_F`$ and $`k_{FIT}=k_{max}`$. At each iteration step, the potential $`U(k)`$ coming directly from the G-matrix calculation is fitted with a parabola, which is then used as the entry potential for the next iteration. Convergence is reached when both potentials remain stable under this procedure. In this procedure one obtains, therefore, two potentials, one calculated from the G-matrix with the parabolic input, and one from the parabolic fit to this potential. Of course, if the the potential coming directly from the G-matrix were indeed parabolic, the two potentials would closely agree. For the choice $`k_{FIT}=2k_F`$ , in Fig. 1 the two potentials at convergence are displayed. In principle, one can calculate the nuclear matter binding energy from both potentials, but the result will be in general slightly different. As one can see, the fully self-consistent potential, obtained without any fitting procedure, as specified above, does not coincide with anyone of the two previous potentials, and these differences give a quantitative indication of the uncertainty introduced by the parabolic approximation. The corresponding saturation curves are reported in Fig. 2a. The parabolic potential produces a saturation curve in fair agreement with the one reported e.g. in ref. . Around saturation the parabolic approximation introduces a shift in the binding of 1-2 MeV. The two choices for the potentials, discussed above, give different binding, since $`U(k)`$ is not really parabolic, and the fitting procedure introduces necessarily an approximation. Another uncertainty is coming from the choice of $`k_{FIT}`$, as can be seen in Fig. 2b, where the results for $`k_{FIT}=k_{max}`$ are reported. In this case the discrepancy are larger for lower density, since then the potential $`U(k)`$ becomes indeed flatter at momenta below $`k_F`$. A more complete account of the dependence on the fitting range is reported in Fig. 3, where the binding at $`k_F=1.4fm^1`$ is reported as a function of $`k_{FIT}`$.

In all cases the saturation curves appear distorted, and the saturation point shifted. Even if in some cases the saturation point seems to be “improved”, this does not have any physical meaning, since, anyhow, it is mainly a spurious effect.

Completely similar results are obtained with the “old” potential Argonne v14 .

In conclusion, we have shown that the parabolic approximation for the single particle potential $`U(k)`$ in the self-consistent Brueckner scheme introduces an uncertainty of 1-2 MeV near the saturation density, and therefore it cannot be used in accurate calculations. The full momentum dependence has to be retained, which prevents the use of a constant effective mass approximation. However, the uncertainty is not dramatic, and for approximate estimates of nuclear binding it can be useful.

Figure captions

Fig. 1.- Single particle potential as function of momentum. The full circles indicate the results of the fully self-consistent calculation, where the potential is taken at each iteration step as calculated from the Brueckner G-matrix. The solid line is the result of the parabolic approximation. The parabolic potential, used as input for the G-matrix, produces the potential indicated by the squares.

Fig. 2.- Saturation curves as obtained for the fully self-consistent calculation (full circles), from the parabolic fit (triangles) and from the single particle potential obtained from the parabolic input (squares). Panel $`a`$ corresponds to $`k_{FIT}=2k_F`$, panel $`b`$ to $`k_{FIT}=k_{max}`$.

Fig. 3.- Dependence of the binding energy from the value of $`k_{FIT}`$, using either the parabolic fit potential (triangles), or the parabolic input (squares) potential. The arrows indicate the value of the binding energy obtained in the fully self-consistent calculation, where no fitting procedure is introduced.


text:
# Creating Spectral Templates from Multicolor Redshift SurveysBased on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

## 1 Introduction

With the introduction of large-format CCDs, the study of galaxy evolution is undergoing a renaissance. Statistically significant volumes of the Universe can now be surveyed at high and low redshift with relatively small amounts of observing time. In fact the rate at which we can undertake deep, multicolor photometric surveys is now over 2 orders of magnitude larger than the rate at which we can followup these observations spectroscopically. As an example, the Sloan Digital Sky Survey will spend 85% of its five year mission measuring the redshifts of one million objects with a unique, multi-fiber spectrograph, while in the remaining 15% of the time it will gather 5 color photometry for 100 million objects.

This rapid growth in multicolor surveys has led to the development of techniques for exploiting the information contained within the photometric data themselves (without relying on followup spectroscopy). One of the most successful of these techniques has been the estimation of galaxy redshifts directly from multicolor data (photometric redshifts) (Baum, 1962). (For an overview of recent work in this field, see Weymann et al. (1999).) These techniques can be broken down into two basic classes. One class of methods use a training set of galaxies with known photometry and spectroscopic redshifts to derive an empirical relation between the colors of galaxies and their redshifts (Connolly et al., 1995a; Wang, Bahcall & Turner, 1998). After calibration, the redshift estimator can predict redshifts for objects with only photometric observations. The advantage of this method lies in its empirical nature; all the effects of dust and galaxy evolution that are present in the training set are implicit within the derived correlation. The drawback is that it can only robustly predict redshifts for objects that are similar to those in the original training set. In other words the estimator is good for interpolation, but not for extrapolation to galaxies that have redshifts or spectral types far from the range covered by the training set.

The second class of techniques, template based photometric redshift estimators (or spectral energy distribution fitting), use model or empirical galaxy spectral energy distributions (Koo, 1985; Connolly et al., 1995b; Gwyn & Hartwick, 1996; Fernández-Soto et al., 1999; Benítez, 1999). By minimizing the difference between the observed and synthetic colors one can find the most likely type and redshift of the galaxy. This method is simple to implement, does not require a training set, and has no extrapolation problems. Its main limitation is that the underlying spectral energy distributions (SEDs) of galaxies within the sample must be well known. Comparisons of the colors of galaxies with spectral synthesis models (Bruzual & Charlot, 1993) have shown that the modeling of the ultraviolet spectral interval for galaxies is uncertain. Whether this is due to uncertainties in the modeling of the stars or due to the effects of dust is as yet unclear. Consequently, photometric-redshift estimates are most accurate when we apply empirical spectral energy distributions derived from observations of local galaxies, e.g. Sawicki et al. (1997). These empirical relations are, however, constructed from a small number of local galaxies that have been observed in detail and there is no guarantee that they represent the full distribution of galaxy types (particularly when we include the effects of evolution with redshift).

In this paper we derive a method that combines the advantages of the above techniques without their drawbacks. We use a training set, not for deriving a direct relationship between colors and redshifts but instead to build an optimal set of spectral templates. We optimize for the shape of the spectral templates to give the best match between the predicted galaxy colors (calculated using the spectroscopic redshifts) and the observed colors. As a result we derive a set of spectral templates that are a better match to the SEDs of the galaxies in the training set than the initial model/empirical templates. Our tests show that these improved templates give a tighter photometric redshift relation than do the original SEDs.

In section 2 we revisit different photometric redshift estimator techniques and we address the question of how to recover information beyond just redshift from photometric observations and how to interpret these results. We introduce an interpolation scheme for optimal spectral type determination. In section 3 we describe our novel technique for creating eigenspectra in an iterative way, and in section 4 we show a test application to the HDF-N/NICMOS catalog, where there are around 150 spectroscopic redshifts available. We present the comparison of photometric redshift estimates based on the most commonly used Coleman, Wu & Weedman (1980) spectra (hereafter CWW) and our empirically developed templates.

## 2 Photometric redshifts

Estimating the redshift of a galaxy from broadband photometry relies on a comparison between the colors of the galaxy and those predicted from either a set of spectral energy distributions (SEDs) or an empirically derived relation between color and redshift. The accuracy of the resulting relation relies critically on having a well calibrated relation (e.g. the spectral templates must cover the full range of galaxy types). In this paper we describe a new technique for developing a set of spectral templates that are optimized to the data in hand. Initially, we outline a variant on the standard template fitting technique that uses an orthogonal basis to compactly describe the broad range of galaxy spectral types.

### 2.1 Data

Unique photometric measurements in four bands (U300, B450, V606, I814) have been available from the Hubble Space Telescope’s Wide Field Planetary Camera 2 (WFPC2) since 1996, when a slice of the very distant Universe became visible in the Hubble Deep Field (HDF-N)(Williams et al., 1996). This data set has been the testbed of photometric redshift estimators from the beginning. Recently the HDF-N field gained two more high signal-to-noise infrared bands (J110, H160) when the Near Infrared Camera and Multi Object Spectrometer (NICMOS) targeted the same area (Dickinson et al., 2000). The joint photometric catalog of the observations was extracted by the SExtractor package (Bertin & Arnouts, 1996) using a weighted sum of the J and H NICMOS images for object detection. In the larger view of the NICMOS instrument, 1681 objects were detected, and out of these, 1291 objects were extracted with complete photometric information in all six bands. In addition to the HST data, we have used $`K_s`$–band photometry derived from images taken at the KPNO 4m telescope (Dickinson, 1998). We have used a method similar to that of Fernández-Soto et al. (1999) for reliably extracting fluxes and errors for every object detected in the NICMOS+WFPC2 images (Papovich et al., 2000). This ensures photometric completeness at $`K_s`$ and eliminates uncertainties concerning aperture corrections between the HST and ground–based data sets. Detection in both optical and infrared enables photometric redshift estimators to identify very high redshift objects out to about $`z=1314`$, because of the continuum blanketing from Lyman-$`\alpha `$ forest (Madau, 1995).

The number of spectroscopic redshifts in HDF-N/NICMOS field is also very promising for various astrophysical projects. Currently there are about 150 spectroscopic redshifts available for the central HDF (Cohen et al., 1996; Steidel et al., 1996; Lowenthal et al., 1997; Dickinson, 1998; Hogg et al., 1998; Weymann et al., 1998; Spinrad et al., 1998; Stern & Spinrad, 1999). These spectroscopic measurements can be used as a training set in empirical photometric redshift estimations or as a test of the accuracy and reliability of template fitting algorithms.

### 2.2 Using eigenspectra

Traditional SED fitting methods provide a simple mechanism for estimating the redshift and spectral type (or metallicity, age etc) of a galaxy. One can compare the measured and the template based colors as a function of redshift and select the redshift and the template where the match is the best. The limitation of this method is that it does not allow a true estimate of the uncertainties on these measures of age etc., nor does it facilitate the propagation of these errors to subsequent analyses of the galaxy distributions. The reason for this is that the multidimensional space of parameters that describe a galaxy are only sampled over a few discrete points (i.e. the number of template spectra).

Instead of using a discrete set of SED templates one can create from them a continuous manifold. We can imagine SED templates as vectors in a multi-dimensional vector space; the dimension of this space ($`D`$) is the number of the wavelength bins of the SEDs. With an orthogonal basis within this vector space one can expand any SED, $`S(\lambda )`$, in terms of this basis,

$$S(\lambda )=\underset{j=1}{\overset{D}{}}c_j\mathrm{\Psi }^j(\lambda )$$
(1)
where $`\mathrm{\Psi }^j(\lambda )`$ is the $`j`$th basis spectrum and $`c_j`$ are the linear expansion coefficients.

One ends up with $`D`$ coefficients to define a SED instead of a single discrete spectral type parameter that selects the template in the traditional way. In practice, the SEDs do not have arbitrary shapes: they occupy only a small a subspace of the possible $`D`$-dimensional vector space. With the help of the the Karhunen-Loève (KL hereafter) transformation (Karhunen, 1947; Loève, 1948; Connolly et al., 1995b; Csabai et al., 2000), a variant of the principal component analysis, we can systematically define a small number of basis vectors ($`\mathrm{\Psi }^j`$), which can reproduce all real spectra to a high (and quantifiable) degree of accuracy. Using only $`\eta `$ eigentemplates the approximated spectrum $`\widehat{S}(\lambda )`$ can be expressed in the form of

$$\widehat{S}(\lambda )=\underset{j=1}{\overset{\eta }{}}c_j\mathrm{\Psi }^j(\lambda )$$
(2)
where the value of $`\eta `$ is typically $`2`$ to $`4`$. We call this an eigenbasis, and its vectors are the eigenspectra. In this way the number of free parameters for which we have to optimize remains small, but we have the freedom of having a continuous distribution of template spectra. Photometric redshift estimation follows the same approach of comparing the derived and actual colors but the coefficients of the eigenspectra are saved as well as the redshift.

Truncating the original basis introduces an uncertainty as a function of wavelength into the linear combination (see equation 2), which can expressed in terms of the remaining eigenspectra as

$$\sigma ^2(\lambda )=\underset{j=\eta +1}{\overset{D}{}}\mathrm{\Lambda }_j\left(\mathrm{\Psi }^j(\lambda )\right)^2$$
(3)
where $`\mathrm{\Lambda }_j`$ is the eigenvalue of the $`j`$th eigenspectra.

Traditionally the eigentemplates are derived from empirical (Coleman, Wu & Weedman, 1980; Kinney et al., 1996) or synthetic (Bruzual & Charlot, 1993) spectra. Unfortunately both methods have limitations when it comes to applications to real data. Empirical spectra are typically derived from observations of a small number of local galaxies. The derived basis may not, therefore, be tenable for high redshift galaxies. A second problem is that the there is only a small range in wavelength where good signal-to-noise data can be obtained; e.g. UV observations or IR spectrophotometry of local galaxies are extremely limited. In contrast, synthetic spectra are available over a wide wavelength range and can be derived for galaxies with a wide variety of ages and past star formation histories, yet many previous studies, e.g. Sawicki et al. (1997); Yee (1998); Benítez (1999), have found that they do not yield photometric redshift estimates as reliable as those derived from empirical templates. Reliable templates could greatly enhance the range where photometric redshift estimation techniques can be applied with acceptable errors.

The eigenbasis approach is more sophisticated in several ways than the traditional discrete template fitting version. As described previously, a small number of continuous parameters are used to determine the rest-frame spectrum, allowing proper error propagation. These parameters are the linear coefficients of the eigenspectra and their combination can be related not only to a single spectral type parameter (as in the traditional version) but also, in principle, to other physical parameters such as reddening, age of the galaxy etc.

While there remains no standard mechanism for defining the spectral type of a galaxy from the eigentemplates, it has been shown that two or three eigentemplates can yield reasonably accurate approximations to most normal galaxy spectra (Connolly et al., 1995b). Both the luminosity of a galaxy and its spectral type are encoded within these coefficients. In this study we use three eigenspectra. Extracting the apparent luminosity from the linear coefficients is straight forward (e.g. the length of the coefficient vector). The normalized Cartesian coordinates can then be converted to a polar system where the spectral type is defined by two mixing angles. The solid line in figure 1 shows these mixing angles as defined from the CWW spectral energy distributions. The galaxy spectral types are found to lie on a one parameter sequence within this two dimensional space.

The points in figure 1 show the mixing angles as defined from the HDF-N WFPC2 and NICMOS data. The ellipses are indicative of the errors associated with the derived mixing angles. Many of the galaxies within the HDF-N data set are seen to follow the locus defined by the CWW templates. There are however a number of outlier data points (with correspondingly large error ellipses). These outliers are typically high redshift objects with reasonably accurate predicted redshifts but with poorly estimated spectra. What happened? The rest-frame spectrum of a high redshift object is only sampled at short wavelengths, thus the cost function in the optimization measuring the deviation of the actual and template based colors is not sensitive to the longer wavelength end. In other words the cost function is “flat” in the expansion coefficient subspace, so tiny deviations drive the resulting spectral type, while the cost vs. redshift has a normal minimum. This observation enables us to introduce a method to avoid negative spectra, while also reducing the number of fitted parameters.

### 2.3 Optimal subspace filtering

We propose a physical interpolation scheme to resolve the type degeneracy of high redshift objects and to insure plausible spectra for any galaxy. We derive a one-dimensional type parameter. As a starting point we use the empirical CWW spectra, extrapolated in the ultraviolet and infrared using empirical templates from Kinney et al. (1996) and models from Bruzual & Charlot (1993), respectively, and with the addition of a model SED (Bruzual & Charlot, 1993) representing a very blue star forming galaxy. We then use these to derive our initial eigenbasis, truncated to a 3-dimensional subspace. The KL expansion of a spectrum on this eigenbasis gives coefficients that represent a galaxy spectrum. We found that the normalized coefficients of the CWW spectra are in a region close to a small circle on the unit sphere. This region defines a one-dimensional sequence, that can be mapped to the $`[0,1]`$ interval, where the value $`0`$ roughly corresponds to “Ellipticals” and $`1`$ represents the bluest galaxies. Deriving the above type parameter consists of the following steps. In the 3D coefficient space, fitting a plane to the normalized coefficient vectors computed for the reference (CWW) spectra gives the equation of a cone, which can be parameterized with one variable. This periodic variable is restricted to the physically sensible interval where the non-negative spectra are. Choosing “0” for the best fit ellipticals and “1” for the bluest galaxies yields an easy interpretation of the estimated type.

Exploiting the above definition of a scalar type parameter, we develop an improved redshift and type estimator (Connolly et al., 1999). In this application there are three continuous fitting parameters: redshift, type and luminosity. In figure 1 open circles represent the mixing angles of the CWW spectra and the solid line shows the projected trajectory from one CWW type to the next. The points in figure 1 show the mixing angles for the HDF sample and the ellipses illustrate the corresponding errors on the derived angles. The type and redshift errors are computed from the shape of the cost function around the minimum, and the formal covariance matrix is obtained from the coefficients of a fitted paraboloid. Looking at figure 1, one can see that the error ellipses on the mixing angles are larger as we move farther from the CWW locus and that the position angles of these error ellipses point in a direction that is almost orthogonal to the CWW curve.

The effect of the optimal subspace filtering is to force the mixing angles to reside on the CWW trajectory. In figure 1 thin lines connect sample objects with error ellipses to the corresponding constrained mixing angles lying on the CWW locus. This physically motivated model—smoothly interpolating between the CWW SEDs—provides sensible spectral types for higher redshift galaxies by suppressing the effect of noise on the derived expansion coefficients. Simply using only the discrete set of four CWW templates would clearly result in unrealistic estimated types for most of the galaxies, eventuating false photometric redshifts.

## 3 Template reconstruction algorithm

We have seen previously that a small number of basis spectra can be used to approximate any spectrum reasonably well. The limitation at this point is the quality of the libraries that the eigentemplates are built from. In order to improve the accuracy of the photometric redshifts and spectral typing we need to build better eigenspectra. The first step was Csabai et al. (2000) using a direct technique to optimize for the shape of the eigenspectra simultaneously with the coefficients. Here we follow a different statistical approach allowing a higher wavelength resolution even for the currently available small data sets (Budavári et al., 1999).

### 3.1 The idea: spectra from broadband photometry

How can we obtain spectra from broadband photometry? We observe galaxies at different redshifts, which means the filters sample the rest-frame light at different wavelengths. If we have large multicolor redshift surveys, then the wavelength regions sampled by the filters (in the rest-frame) overlap. Figure 2 illustrates how the HDF/NICMOS photometric bands (projected back to the rest-frame) overlap for a set of spectroscopic redshifts taken from the catalog. Looking at a particular wavelength, we can actually see which bands set up constraints on the value of the eigenspectra. This oversampling of the eigenspectra allows us to develop applications to extract higher resolution information. Developing such a reconstruction process is not trivial, though. Several problems emerge because of the non-uniform redshift distribution and the photometric uncertainties within the data. Obviously none of the standard de-convolution applications can be used to solve the problem. We have, therefore, developed an iterative procedure for correcting eigenspectra, which repairs the eigenspectra only where sufficient information exists and where the magnitude of the modification is determined statistically.

Our iterative procedure starts from an initial set of eigentemplates and creates increasingly better eigenbases by improving the match between the measured colors and those derived from the templates. The algorithm is quite simple. Given a catalog of galaxies with broadband photometry and spectroscopic redshifts, we first expand each galaxy over our initial set of templates, transformed to the galaxy’s redshift, by solving for the linear expansion coefficients.

If the initial templates do not accurately reproduce the ensemble of galaxies the spectrum of every galaxy will have systematic errors due to the choice of templates. This corresponds to a subspace in the space of all possible spectra. A better representation of the photometry may require moving out of this subspace. We can achieve this by an iterative repair of the estimated spectrum of each galaxy, by considering what is the minimal change in its shape that would improve the photometry. We compute this “repair” by minimizing a cost-function over the shape of the spectrum, penalizing both large deviations from the photometry and from the template expansion. Since the template expansion was already an optimal representation given the subspace defined by the templates, this repair enables the spectra to move out of the original subspace.

Building a new eigenbasis over the set of repaired spectra will allow changes on the templates which improve the fit population-wide, but reject individual fluctuations. Performing this process multiple times in an iterative fashion (Everson & Sirovich, 1995), the overall fit to the photometry improves substantially. The flowchart of the algorithm is shown in figure 3.

### 3.2 Repairing spectra

The key step to the template reconstruction algorithm is the spectrum repairing method. The spectrum we are about to repair is an approximation of the real spectrum based on the photometric data. This best fitting spectrum can be found within the subspace spanned by the current eigenspectra by varying the coefficients directly. To obtain a physically sensible spectrum the optimal subspace filtering method (see section 2.3) should be used. Having a set of eigentemplates and the calculated best fitting type the corresponding rest-frame spectrum of a galaxy can be easily computed. The best fitting type $`t`$ is calculated by minimizing the cost function,

$$C^2(t)=\underset{n}{}\frac{1}{\mathrm{\Delta }_n^2}\left(\widehat{f}_nf_n(t,z_{\mathrm{spec}})\right)^2$$
(4)
where $`\widehat{f}_n`$, $`\mathrm{\Delta }_n`$ are the measured flux and error in the $`n`$th passband. In fact, the type parameter $`t`$ determines the values of the coefficients of the eigenspectra and the linear combination gives the best fitting spectrum. Using $`\eta `$ number of eigentemplates the approximated spectrum can be expressed in the form of

$$\widehat{S}(\lambda )=\underset{j=1}{\overset{\eta }{}}c_j(t)\mathrm{\Psi }^j(\lambda )$$
(5)
where $`\mathrm{\Psi }^j(\lambda )`$ is the $`j`$th eigenspectrum and $`c_j(t)`$ are the linear coefficients defined by the type $`t`$ and the brightness.

The basic idea is to modify this spectrum to match the observed colors better. This is another minimization problem, where we try to adjust the spectrum based on a cost function built up by two different terms. On one hand, we have the deviation of the actual photometry from the fluxes computed from the spectrum itself, $`S(\lambda )`$. On the other hand we have the deviation of the spectrum $`S(\lambda )`$ from the template expansion $`\widehat{S}(\lambda )`$. The repair process has to be able to deal properly with different types of errors such as the photometric errors and the uncertainties within the eigenspectra.

Let $`R^n(\lambda )`$ represent the response function of the $`n`$th filter, normalized and corrected for the detector and telescope throughput function. The derived flux in the $`n`$th passband can be written as an integral or approximated with the sum over the discrete representation of the functions.

$$f_n=R^n(\lambda )S(\lambda )𝑑\lambda =\underset{k}{}r_k^ns_k$$
(6)
where $`s_k`$ and $`r_k^n`$ are the discrete representation of the spectrum $`S(\lambda )`$ and the filter $`R^n(\lambda )`$ and $`k`$ refers to a certain wavelength, $`\lambda _k`$. The discrete version of $`\widehat{S}(\lambda )`$ is given by the linear combination of $`c_j(t)`$ and $`\psi _k^j`$, the new representation of the $`j`$th eigenspectrum $`\mathrm{\Psi }^j(\lambda )`$.

The discrete problem translates into a multidimensional minimization of the cost function as a function of the unknown spectral shape $`s_k`$.

$$\chi ^2=\underset{k}{}\frac{1}{\sigma _k^2}\left(s_k\widehat{s}_k\right)^2+\underset{n}{}\frac{1}{\mathrm{\Delta }_n^2}\left(f_n\widehat{f}_n\right)^2$$
(7)
where $`\sigma _k`$ describes the ability of the eigenspectra to be changed at a certain wavelength. The primary source of error in the expansion is the truncation error due to the restricted subspace; see equation 3 and its discrete form,

$$\sigma _k^2=\underset{j>\eta }{}\mathrm{\Lambda }_j\left(\psi _k^j\right)^2.$$
(8)
The cost function of the problem can be transformed into a simpler form by introducing the difference spectrum ($`\xi `$) and some constant terms ($`g_n`$).

$`\xi _k`$ $`=`$ $`s_k\widehat{s}_k`$ (9)
$`g_n`$ $`=`$ $`\widehat{f}_n{\displaystyle \underset{k}{}}r_k^n\widehat{s}_k`$ (10)
$$\chi ^2=\underset{k}{}\frac{\xi _k^2}{\sigma _k^2}+\underset{n}{}\frac{1}{\mathrm{\Delta }_n^2}\left(g_n\underset{k}{}r_k^n\xi _k\right)^2$$
(11)
The cost function is quadratic in its variables, $`\{\xi _k\}`$. The minimum can be analytically derived, since the partial derivatives vanish at the minimum. This yields a system of linear equations, which can be solved for $`\xi `$.

$$\frac{\chi ^2}{\xi _l}=0;\underset{k}{}M_{lk}\xi _k=\nu _l$$
(12)
where

$$M_{lk}=\underset{n}{}\frac{r_l^nr_k^n}{\mathrm{\Delta }_n^2}+\frac{\delta _{lk}}{\sigma _k^2}\mathrm{and}\nu _l=\underset{n}{}\frac{r_l^ng_n}{\mathrm{\Delta }_n^2}$$
(13)
and $`\delta _{lk}`$ is the usual Kronecker symbol. The repaired spectrum can be computed from the difference vector by simply adding it to the template based spectrum, $`s_k=\widehat{s}_k+\xi _k`$. A few typical example repaired spectra can be seen in figure 4 for objects with different spectral type parameters. The steps we take are small because the role of these tiny deviations from the template based spectra is to point the eigenspectra into the right direction in a statistical sense.

Once the spectrum repairing procedure is applied to each galaxy in the training catalog, the Karhunen-Loève transformation can be applied to obtain a new set of eigenspectra. Repeating the above steps yields a statistically robust training algorithm to develop empirical eigenspectra from photometry. If there were no modifications to the linear combinations of the eigenspectra then we would get exactly the same spectrum subspace at the end. Simply adding some uncorrelated noise to the spectra would retain the subspace as well. Corrections in the eigenspectra can only occur due to changes correlated with galaxy type.

## 4 Application to the HDF/NICMOS data

### 4.1 Building the templates

Our algorithm is constructed in such a way that the iteration can be initialized with any eigentemplates and started from imported spectra. In our study for deriving the initial set of eigentemplates we used the extended CWW spectra, which emerged as the most commonly used set of templates by several other groups.

The actual implementation was somewhat different. Logarithmically resampled templates were used in order to speed up the computation, since a logarithmically rebinned spectrum can be redshifted by simply offsetting the indices of the vectors. The result of this is that the integral measure changes, which requires some trivial modifications in the above equations. Having more bins in the blue than the red ranges (up to 25000Å) allows higher resolution in the optical bands.

Repairing the eigentemplates improved the rest-frame colors rapidly. Global features changed significantly over the initial iterations. After just a few steps changes to the continuum of the eigenspectra converge. Figure 5 illustrates the continuum changes of the basis after just two iterations. Subsequent iterations change the higher frequency components within the eigenspectra. Figure 6 and table 1 containing the actual numeric values of the estimates help us to understand the effects of these iterations on the estimated spectra. Changes in the eigenbasis result in a new fit, the expansion coefficients differ and thus the actual estimated SEDs can differ substantially more than one might expect from figure 5. Consistent with our expectations about the uncertainties present within the empirical SEDs, most modifications occur in the ultraviolet and the infra red, where the initial templates are poorly constrained. The bluest starburst SEDs tend to be changed least.

Each iteration requires approximately 2 minutes on a workstation with the CPU time scaling linearly with the number of objects (as most of the time was spent in the spectrum correction phase solving the linear equations for the modified spectrum). The KL expansion took an insignificant amount of time on such a small data set. However, even the principal component analysis can be done iteratively, which enables us to apply this technique to much bigger surveys, such as the Sloan Digital Sky Survey. Our training algorithm is applicable to data sets composed of catalogs in different photometric systems. Theoretically there is nothing to prevent us from including all observed galaxies which have photometry and redshifts available.

The stability of our technique had been tested with randomly chosen subsets of the entire training set. The iterated eigenbases were compared to the original and the difference of their projections onto the same subspace has been analyzed. Calculating the mean and its standard deviation of thirty resulting templates in all the wavelength bins shows us that these changes are significant. Figure 7 show the mean difference of the resulting and CWW eigentemplates with solid lines and the scatter with dotted lines. The first eigencomponent changes less but it is the most significant.

### 4.2 Applying to photometric redshifts

In our analysis the spectroscopic redshifts have only been used to shift the templates to a common rest-frame for the training set. The template repairing algorithm does not yield photometric redshifts nor is it directly related to the estimation at all. It reconstructs the (continuum) rest-frame shape of galaxy SEDs. These eigenspectra can be utilized in many different ways; to use them as templates in photometric redshift estimation is just one possibility.

Figure 8 shows how the redshift prediction improves as the eigenspectra evolve. The estimates are approaching the ideal $`z_{\mathrm{phot}}=z_{\mathrm{spec}}`$ line as the SEDs converge to the best possible representation of reality. Figure 8 is a comparison of the photometric redshift prediction based on the resulting KL (bottom) and CWW (top) eigentemplates. The systematic errors disappear and the scatter reduces. Beyond the standard root mean square error ($`\mathrm{\Delta }_{\mathrm{rms}}`$) the relative deviation of $`(1+z)`$ was also computed ($`\mathrm{\Delta }_{\mathrm{rel}}`$) for comparison with estimates found by other authors.

$$\mathrm{\Delta }_{\mathrm{rms}}=\left(z_{\mathrm{spec}}z_{\mathrm{phot}}\right)^2$$
(14)
$$\mathrm{\Delta }_{\mathrm{rel}}=\left(\frac{z_{\mathrm{spec}}z_{\mathrm{phot}}}{1+z_{\mathrm{spec}}}\right)^2$$
(15)
Table 2 compares the overall redshift errors computed for both the CWW and the trained (KL) eigenspectra. It shows significant improvements on any redshift range. The final KL eigenbasis gains about a factor of two in the average regardless of which statistics are preferred. The rms error changes from 0.23 to 0.12 for estimates up to redshift of 6.

An important feature of this approach is the rather rigorous handling of error propagation. The available photometric errors were incorporated in both the template repair and photometric redshift estimation routines. The photometric inversion should predict not only the physical parameters but also their errors (that are usually correlated). Since it is hard to find the global minimum of a function in more dimensions, we evaluate the cost function on a grid to speed up the processing and to ensure the optimization does not end up in a local minimum. Searching the grid on the type-redshift plane provides a good estimate of the global minimum, however the results can be revised using the continuous cost function again. Fitting a second order surface around the minimum allows us to refine the estimation by choosing the fitted minimum and to compute the formal covariance matrix from the coefficients of the fitting formula. Figure 9 shows the correlated type and redshift errors for the training set. Consistent with expectations, most of the ellipses show that higher redshifts than those estimated are possible if we choose an incorrect (bluer) template. This illustrates the need to fully and finely sample the distribution of template spectra. The bluer SEDs also tend to have larger errors on the type and redshift than the ellipticals due to the less significant continuum features in their spectra.

## 5 Discussion

Our technique has been shown to be able to derive rest-frame templates from photometry in an iterative way, improving step-by-step the correlation of the template based colors to the actual measurements. The statistical approach of the KL expansion provides a robust determination of the systematic changes in the SEDs. After just the initial 2-3 iterations, the overall continuum shapes of the SEDs are adjusted to reproduce the photometry. Higher resolution detailed features emerge in later iterations. Consistent with our expectations, significant continuum changes occur both in the ultraviolet and infrared, where the original templates are poorly constrained. The estimated SEDs based on the trained eigenbases tend to be redder in general. The UV part of the spectra get fainter and the flux increases in the IR for the early type galaxies (which change the most). The blue galaxies have gone through less dramatic modulation. Clearly, the color spread of galaxies is larger than one might anticipate from the CWW spectra.

The limitation on our current application is the relatively small number of objects with accurate multicolor photometry and spectroscopic redshifts in the Hubble Deep Field. Combining several observations into one training set is one way to improve the robustness of our results. Our technique is constructed to be able to deal with many multicolor observations at the same time. Data of various surveys can be incorporated regardless of photometric systems they use, meaning not just broad but also narrow or intermediate band imaging. This flexibility of the algorithm allows us to simply merge data sets coming from several instruments, where the data set up constraints on the SEDs at different wavelengths. With time we will be able to improve other parts of the SEDs. New upcoming instruments will provide high signal-to-noise ratio data in different spectral ranges, e.g. the Galaxy Evolution Explorer (GALEX) probing the far and near ultraviolet (1300-3000Å) or the Space Infrared Telescope Facility (SIRTF) and the Next Generation Space Telescope (NGST) looking at the infrared ranges. In fact, the algorithm is even more generic. Eventually, spectrophotometric observations can be naturally incorporated in the training mechanism.

There remains one class of galaxy for which this technique is not optimal; high redshift galaxies for which we cannot easily obtain spectroscopic redshifts can not be optimized by our procedure (i.e. we lack the relevant training set). As our technique can seamlessly incorporate theoretical and empirical data (weighting the samples appropriately) if we have a reasonable understanding of the stellar populations we can account for missing galaxy populations. Objects of extremely unusual colors will not be well matched by any template fitting procedure but the goodness of fit to these data will enable us to identify these anomalous objects.

In this paper we have presented a new method to reconstruct the rest-frame continuum spectra of galaxies directly from a catalog of broadband multicolor photometric observations and spectroscopic redshifts. In addition to an improved photometric redshift estimation the application of eigenspectra to this problem enables the classification of galaxy spectra to move from a series of discrete (and somewhat ad hoc) classification types to a continuous parameterization. With the help of the optimal subspace filtering, non-physical spectral energy distributions could be avoided. As a result of the repair algorithm, the eigenbasis improved the capability to represent the observed photometric properties. This improvement cannot be tested directly, since we do not know the real, photometrically calibrated spectra for the objects in our training catalog, but the photometric redshift estimation test with the repaired spectra gave substantially better results than the estimation using model/empirical templates ($`\mathrm{\Delta }_{\mathrm{rms}}=0.12`$ vs. $`\mathrm{\Delta }_{\mathrm{rms}}=0.23`$). New photometric and spectroscopic surveys will extend the size of the training catalog, making possible to improve further the spectral templates and enrich them with finer details.

IC and TB acknowledges partial support from the MTA-NSF grant no. 124 and the Hungarian National Scientific Research Foundation (OTKA) grant no. T030836, AJC acknowledges support from an LTSA grant (NAG57934) and from NASA through grant number GO-07817.06-96A from the Space Telescope Science Institute. AS acknowledges support from NSF (AST9802980) and a NASA LTSA (NAG53503). MD received support from NASA through grant number GO-07817.01-96A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.