format: text

text:
# Gravitational Wave Damping of Neutron Star Wobble

## I Introduction

This paper calculates the effect of gravitational wave (gw) back-reaction on the torque-free precession, or wobble, of realistic, spinning neutron stars (NS’s). By ‘realistic’ we mean the NS is treated as a mostly-fluid body with an elastic crust, as opposed to a rigid body. (However we do not include any superfluid effects in our analysis.) Freely precessing neutron stars are a possible source for the laser interferometer gw detectors (LIGO, VIRGO and GEO under construction, TAMA already operational); it is the prospect of gravitational wave astronomy that motivated our study. Also, the first clear observation of free precession in a pulsar signal was reported very recently , making this investigation all the more timely.

The effect of gw back-reaction on wobbling, axisymmetric rigid bodies was first derived 27 years ago, in an impressively early calculation by Bertotti and Anile . They found (correctly) that for rigid bodies, gw backreaction damps wobble on a timescale (for small wobble angle $`\theta `$) $`\tau _\theta ^{\mathrm{rigid}}=1.8\times 10^6\mathrm{yr}[10^7/(\mathrm{\Delta }I/I_1)]^2(\mathrm{kHz}/\nu _s)^4(10^{45}\mathrm{g}\mathrm{cm}^2/I_1`$), where $`\nu _s`$ is the spin frequency and $`\mathrm{\Delta }I=(I_3I_1)`$ (with $`I_1=I_2I_3)`$.

In the same paper, Bertotti and Anile  went on to calculate the effect of gw back-reaction on wobble for the more realistic case of an elastic NS. When cast into our notation, their claimed gw timescale is $`5I_1c^5/[2G(2\pi \nu _s)^4\mathrm{\Delta }I_\mathrm{\Omega }\mathrm{\Delta }I_d]`$ where $`\mathrm{\Delta }I_\mathrm{\Omega }`$ is the asymmetry in the moment of inertia due to centrifugal forces and $`\mathrm{\Delta }I_d`$ is the asymmetry due to some other mechanism, such as strain in the solid crust. Taking $`\mathrm{\Delta }I_\mathrm{\Omega }`$ to be (roughly) the asymmetry expected for a rotating fluid according to $`\mathrm{\Delta }I_\mathrm{\Omega }/I0.3(\nu _s/\mathrm{kHz})^2`$, we would then have a damping time of merely $`0.6\mathrm{yr}(\mathrm{kHz}/\nu _s)^6[10^7/(\mathrm{\Delta }I_d/I)](10^{45}\mathrm{g}\mathrm{cm}^2/I)`$. Despite the fundamental beauty of this probem and its potential astrophysical significance, their remarkable claim–that in realistic NS’s, gw’s damp wobble with amazing efficiency–was apparently little known. (A citation index search showed that Bertotti and Anile  had been referenced by other authors only four times in the last 27 years.)

We will show that the Bertotti and Anile result for elastic NS’s is very wrong, however. For typical cases of interest, their gw timescale $`\tau _\theta `$ is too short by a factor $`10^5`$ ! Moreover, their calculation even gives the wrong sign (exponential growth instead of damping) when $`\mathrm{\Delta }I_d`$ is negative. ***Actually, Bertotti and Anile  never claim in words that they find unstable growth of the wobble angle when $`\mathrm{\Delta }I_d<0`$, but that is what is found if one just takes their formulae and converts from their notation to ours, as above. Moreover we have repeated their (flawed) calculation, including their one crucial error, and seen that it does lead to a prediction of exponential wobble growth for $`\mathrm{\Delta }I_d`$ negative. The conversion from their notation to ours is simply $`(\delta _1I\delta _2I)(\mathrm{cos}^2\gamma \frac{1}{2}\mathrm{sin}^2\gamma )\mathrm{\Delta }I_d`$ and $`\delta _2I\mathrm{\Delta }I_\mathrm{\Omega }`$). In contrast, we find that gw’s always act to damp the wobble in realistic NS’s, just as for rigid bodies. While in Nature the typical case will be $`\mathrm{\Delta }I_d`$ positive, $`\mathrm{\Delta }I_d<0`$ can also occur in principle. We call attention to this case not because it is common, but because it highlights how much our result differs from Bertotti and Anile , and because, in fact, their implicit prediction of exponential wobble growth for this case provided our initial impetus to look more closely at this problem.

The organization of this paper is as follows. In §2 we derive the gw damping timescale for rigid-body wobble, using the mass quadrupole expressions for the energy and angular momentum radiated to infinity. (This derivation is actually Exercise 16.13 in the textbook by Shapiro and Teukolsky ). We give another derivation of $`\tau _\theta `$ in §3, this time by adding the Burke-Thorne radiation reaction force directly to the Newtonian equations of motion. This latter approach was how Bertotti and Anile  first calculated (correctly) the gw damping time for wobbling, rigid bodies.

In §4 we review standard material on the torque-free precession of elastic bodies, in the absence of viscous terms or gw back-reaction. In §5 we derive the gw damping timescale $`\tau _\theta `$ in the elastic case, using energy and angular momentum balance. In §6 we give a second derivation of $`\tau _\theta `$ in the elastic case, using the Burke-Thorne radiation reaction force to evolve the elastic body’s free precession. This was also the strategy of Bertotti and Anile , and we show where they went wrong. Briefly, they did not realize that in addition to torquing the NS, the radiation reaction force also perturbs the NS’s shape (in particular, its inertia tensor). When solving for the evolution of the wobble angle, we show that the “perturbed shape” term in the equations of motion almost entirely cancels the gw torque term that they do include. (Of course, by definition there is no “perturbed shape” term in the rigid-body case, which is probably why they forgot this term when adapting that calculation to the elastic case.) In §7 we describe how to include the effects of a fluid core in the radiation reaction calculation. Finally, in §8 we conclude by commenting briefly on the astrophysical implications of our result.

We will work in cgs units.

## II Radiation reaction for a rigid body: Energy and angular momentum balance

The derivation of the wobble damping rate for realistic NS’s, using energy and angular momentum balance, is rather similar to the corresponding derivation for rigid bodies. Here we briefly review the solution to the rigid-body problem, as a warm-up for tackling the realistic case.

Consider an axisymmetric rigid body with principal axes $`\widehat{x}_1,\widehat{x}_2,\widehat{x}_3`$ and principal moments of inertia $`I_1=I_2I_3`$. Let the body have angular momentum $`𝑱`$, misaligned from $`\widehat{x}_3`$. Define the wobble angle $`\theta `$ by $`𝑱\widehat{𝐱}_\mathrm{𝟑}=J\mathrm{cos}\theta `$. It is a standard result from classical mechanics that (in the absence of external torques) the body axis $`\widehat{x}_3`$ precesses around $`𝑱`$ with (inertial frame) precession frequency $`\dot{\varphi }=J/I_1`$, with $`\theta `$ constant . Together, the pair $`(\theta ,\dot{\varphi })`$ completely specify the free precession (modulo a trivial constant of integration specifying $`\varphi `$ at $`t=0`$). We wish to calculate the evolution of these two parameters using the *time-averaged* fluxes $`(\dot{E},\dot{J})`$.

Straightforward application of the mass quadrupole formalism  gives

$$\dot{E}=\frac{2G}{5c^5}\dot{\varphi }^6(\mathrm{\Delta }I)^2\mathrm{sin}^2\theta (\mathrm{cos}^2\theta +16\mathrm{sin}^2\theta ),$$
(1)
where $`\mathrm{\Delta }I=I_3I_1`$, and

$$\dot{J}=\dot{E}/\dot{\varphi }.$$
(2)
It follows from differentiation of $`\dot{\varphi }=J/I_1`$ that

$$\ddot{\varphi }=\frac{2G}{5c^5}\frac{\mathrm{\Delta }I^2}{I_1}\dot{\varphi }^5\mathrm{sin}^2\theta (16\mathrm{sin}^2\theta +\mathrm{cos}^2\theta ).$$
(3)
To calculate the rate of change of the wobble angle, rearrange

$$\frac{dE}{dt}=\frac{E}{J}|_\theta \frac{dJ}{dt}+\frac{E}{\theta }|_J\frac{d\theta }{dt}$$
(4)
to give

$$\dot{\theta }=\frac{\dot{J}[\dot{\varphi }\frac{E}{J}|_\theta ]}{\frac{E}{\theta }|_J},$$
(5)
where Eq. (2) has been used. The energy of the body is simply its kinetic energy:

$$E=\frac{J^2}{2I_1}\left[1\mathrm{cos}^2\theta \frac{\mathrm{\Delta }I}{I_3}\right],$$
(6)
and so

$$\frac{E}{J}|_\theta =\frac{J}{I_1}\left[1\mathrm{cos}^2\theta \frac{\mathrm{\Delta }I}{I_3}\right],$$
(7)
$$\frac{E}{\theta }|_J=\frac{J^2}{I_1}\mathrm{cos}\theta \mathrm{sin}\theta \frac{\mathrm{\Delta }I}{I_3}.$$
(8)
This gives

$$\dot{\theta }=\frac{2G}{5c^5}\frac{\mathrm{\Delta }I^2}{I_1}\dot{\varphi }^4\mathrm{cos}\theta \mathrm{sin}\theta (16\mathrm{sin}^2\theta +\mathrm{cos}^2\theta ).$$
(9)
We can construct timescales on which the spin-down and alignment occur:

$$\tau _{\dot{\varphi }}^{\mathrm{rigid}}=\frac{\dot{\varphi }}{\ddot{\varphi }}=\frac{5c^5}{2G}\frac{1}{\dot{\varphi }^4}\frac{I_1}{\mathrm{\Delta }I^2}\frac{1}{\mathrm{sin}^2\theta (16\mathrm{sin}^2\theta +\mathrm{cos}^2\theta )},$$
(10)
$$\tau _\theta ^{\mathrm{rigid}}=\frac{\mathrm{sin}\theta }{\frac{d}{dt}\mathrm{sin}\theta }=\frac{5c^5}{2G}\frac{1}{\dot{\varphi }^4}\frac{I_1}{\mathrm{\Delta }I^2}\frac{1}{\mathrm{cos}^2\theta (16\mathrm{sin}^2\theta +\mathrm{cos}^2\theta )}.$$
(11)
Radiation reaction causes both $`\dot{\varphi }`$ and $`\mathrm{sin}\theta `$ to decrease, regardless of whether the body is oblate or prolate. Note that in the limit of small wobble angle the inertial precession frequency remains almost constant ($`\tau _{\dot{\varphi }}^{\mathrm{rigid}}\mathrm{}`$), while $`\theta `$ decreases exponentially on the timescale

$$\tau _{\theta 1}^{\mathrm{rigid}}=\frac{5c^5}{2G}\frac{1}{\dot{\varphi }^4}\frac{I_1}{\mathrm{\Delta }I^2}.$$
(12)
Parameterising:

$$\tau _\theta ^{\mathrm{rigid}}=1.8\times 10^6\mathrm{yr}\left(\frac{10^7}{\mathrm{\Delta }I/I_1}\right)^2\left(\frac{\mathrm{kHz}}{\nu _s}\right)^4\left(\frac{10^{45}\mathrm{g}\mathrm{cm}^2}{I_1}\right).$$
(13)
In the limit of vanishingly small wobble angle the partial derivative on the lhs of Eq. (7) becomes what we conventionally call the ‘spin frequency’ $`\mathrm{\Omega }`$ of the body . Eq. (5) then shows that $`\dot{\theta }`$ is proportional to the difference between the inertial precession frequency $`\dot{\varphi }`$ and the spin frequency $`\mathrm{\Omega }`$. This difference remains finite as $`\theta 0`$ according to $`\dot{\varphi }\mathrm{\Omega }=(\mathrm{\Delta }I/I_1)\mathrm{\Omega }[1+𝒪(\theta ^2)]`$. Thus for a prolate body ($`\mathrm{\Delta }I<0`$), such as an American football, the body precesses slower than it spins, while for an oblate body the inertial precession frequency is higher than the spin frequency. Since the denominator in (5) is also proportional to $`\mathrm{\Delta }I`$, the wobble angle decreases regardless of the sign of this factor. This viewpoint will be useful when we consider the radiation reaction problem for an elastic body.

## III Radiation reaction for rigid bodies: Local force

We will now re-derive the spin-down and alignment timescales by adding the Burke-Thorne local radiation reaction force to the equations of motion.

The Burke-Thorne radiation reaction potential at a point $`x`$ is given by :

$$\mathrm{\Phi }^{RR}=\frac{G}{c^5}x^ax^b\frac{d^5I\text{ }_{ab}}{dt^5},$$
(14)
where $`I\text{ }_{ab}`$ denotes the trace-reduced quadrupole moment tensor:

$$I\text{ }_{ab}=_V\rho (x_ax_b\frac{1}{3}\delta _{ab}x^2)𝑑V.$$
(15)
Note that this is related to the moment of inertia tensor according to

$$I\text{ }_{ab}=I_{ab}\frac{2}{3}\delta _{ab}_V\rho x^2𝑑V,$$
(16)
with the result that

$$\mathrm{\Delta }II_3I_1=(I\text{ }_3I\text{ }_1).$$
(17)
The radiation reaction force (on a particle of unit mass) is $`F_a^{RR}=\mathrm{\Phi }^{RR}/x^a`$. The instantaneous (*not time-averaged*) torque on a body can easily be shown to be

$$T^a=\frac{2G}{5c^5}ϵ^{abc}I\text{ }_{bd}\frac{d^5I\text{ }_{dc}}{dt^5}.$$
(18)
Making use of Eq. (17) it is straightforward to calculate this torque for the free precessional motion. We find

$$𝑻=\frac{2G}{5c^5}\mathrm{\Delta }I^2\dot{\varphi }^5\mathrm{sin}\theta (16\mathrm{sin}^2\theta +\mathrm{cos}^2\theta )𝒏_{\mathbf{}𝒏_𝒅},$$
(19)
acting always in the plane containing the angular momentum and the symmetry axis $`x_3`$, and perpendicular to $`𝒏_𝒅`$, i.e. along the direction of $`𝒏_{\mathbf{}𝒏_𝒅}`$ shown in Fig. 1. We will refer to this plane as the *reference plane*.

The evolution equations can be calculated without going to the trouble of writing down Euler’s equations. Differentiation of $`\dot{\varphi }=J/I_1`$ gives

$$\ddot{\varphi }=\frac{\dot{J}}{I_1},$$
(20)
and so

$$\ddot{\varphi }=\frac{T\mathrm{sin}\theta }{I_1}.$$
(21)
Define $`J_{n_d}`$ as the component of the angular momentum perpendicular to the symmetry axis. Then differentiation of the trivial relation

$$\mathrm{sin}\theta =\frac{J_{n_d}}{J}$$
(22)
leads to

$$\dot{\theta }=\frac{T_J}{J}=\frac{T\mathrm{cos}\theta }{J},$$
(23)
where $`T_J`$ is the component of the torque perpendicular to $`𝑱`$. Equations (21) and (23) show that the action of the torque breaks down neatly into two parts. The component along $`𝑱`$ acts to change the inertial precession frequency $`\dot{\varphi }`$ while the component perpendicular to $`𝑱`$ acts to change $`\theta `$. Substitution of (19) into (21) and (23) then reproduces the spin-down and alignment of equations (3) and (9), so the two methods of calculation agree. As this torque formulation makes clear (by combining Eqs. 21 and 23), the product $`\dot{\varphi }\mathrm{cos}\theta `$ remains constant, so that if a body is set into free precession described by $`(\theta _0,\dot{\varphi }_0)`$, it tends to a non-precessing motion about $`x_3`$ with (inertial frame) angular velocity $`\dot{\varphi }=\mathrm{cos}\theta _0\dot{\varphi }_0`$.

## IV Torque-free Precession of Elastic Bodies

We now review the theory of the free precession of an elastic body. This problem was first addressed in the context of the Earth’s own motion. A rigorous treatment of the methods employed can be found in Munk and MacDonald . The terrestrial analysis was extended to neutron stars by Pines and Shaham . The energy loss due to gravitational waves was considered by Alpar and Pines .

Following the latter authors we will model a star consisting of a centrifugal bulge and a *single* additional deformation bulge. Alpar and Pines wrote an inertia tensor for the elastic body of the form

$$𝑰=I_{0,S}𝜹+\mathrm{\Delta }I_\mathrm{\Omega }(𝒏_𝛀𝒏_𝛀1/3𝜹)+\mathrm{\Delta }I_d(𝒏_𝒅𝒏_𝒅1/3𝜹).$$
(24)
where $`𝜹`$ is the unit tensor $`[1,1,1]`$, $`𝒏_𝛀`$ is the unit vector along the star’s angular velocity $`𝛀`$, and $`𝒏_𝒅`$ is the unit vector along the body’s principal deformation axis (explained below). The $`I_{0,S}`$ and $`\mathrm{\Delta }I_d`$ pieces of $`𝑰`$ together represent the inertia tensor for the corresponding non-rotating star. The $`\mathrm{\Delta }I_d`$ term is just the non-spherical piece of this tensor (approximated as axisymmetric). If the star were a perfect fluid, $`\mathrm{\Delta }I_d`$ would vanish, but in real stars (and the Earth) $`\mathrm{\Delta }I_d`$ is non-zero due to crustal shear stresses and magnetic fields. The term $`\mathrm{\Delta }I_\mathrm{\Omega }`$ ($`>0`$ and $`\mathrm{\Omega }^2`$ for small $`\mathrm{\Omega }`$), represents the increase in the star’s moment of inertia (compared to the non-rotating case) due to centrifugal forces. Since the crust of a rotating NS will tend to “relax” towards its oblate shape, having $`\mathrm{\Delta }I_d>0`$ is surely the typical case in Nature. (E.g., if one could slow the Earth down to zero angular velocity without cracking its crust, it would remain somewhat oblate: the crust’s “relaxed, zero-strain” shape is oblate, and after centrifugal forces are removed, the stresses that build up in the crust will act to push it back towards that relaxed shape.) But a negative $`\mathrm{\Delta }I_d`$ is also possible in principle. We say the deformation bulge aligned with $`𝒏_𝒅`$ is ‘oblate’ if $`\mathrm{\Delta }I_d>0`$ and ‘prolate’ if $`\mathrm{\Delta }I_d<0`$.

What is a typical magnitude for $`\mathrm{\Delta }I_d`$ in real, spinning NS’s? Let us assume $`\mathrm{\Delta }I_d`$ is due primarily to crustal shear stresses (as opposed to stresses in a hypothetical solid core, extremely strong B-fields, or pinned superfluid vortices). Then for a relaxed crust (i.e., a crust whose reference ellipticity is very close to its actual ellipticity), we have $`\mathrm{\Delta }I_d=b\mathrm{\Delta }I_\mathrm{\Omega }`$, where Alpar and Pines  estimate $`b10^5`$ for a primordial (cold catalyzed) crust. The maximum value for $`\mathrm{\Delta }I_d/I`$ is therefore of order $`10^5`$. The parameter $`b`$ (which arises from inter-nucleon Coulomb forces) scales like the average $`Z^2/A`$ of the crustal nuclei. Since crusts of accreted matter (as in LMXB’s) have smaller-Z nuclei , their $`b`$ factor is correspondingly smaller, by a factor $`23`$. Using $`\mathrm{\Delta }I_\mathrm{\Omega }/I0.3(\nu _s/\mathrm{kHz})^2`$, we would therefore estimate $`\mathrm{\Delta }I_d/I10^7`$ for a NS with a relaxed, accreted crust and $`\nu _s300`$ Hz, while for the Crab one would expect $`\mathrm{\Delta }I_d/I3\times 10^9`$ (again, assuming its crust is almost relaxed). For the freely precessing pulsar reported in Stairs et al. , where the body-frame precession period is $`2\times 10^8`$ times the rotation period, Eq. (38) below (valid for elastic bodies) yields $`\mathrm{\Delta }I_d/I=5\times 10^9`$. For $`b=10^5`$ this corresponds to a reference oblateness of $`5\times 10^4`$. This is consistent with the star’s crust having solidified when it was spinning at about 40 Hz, assuming that neither glitches nor plastic flow have modified its shape since. (When the effects of crust-core coupling are taken into account, giving Eq. (65), this initial frequency reduces to 12 Hz. See Jones  for a review of pulsar free precession observations).

Precession occurs when $`𝒏_𝒅`$ and $`𝒏_𝛀`$ are not aligned. Below we describe the precessional motion when there is no damping. This analysis is quite general: it applies to any star whose inertia tensor is described by Eq. (24), independent of what causes the deformation bulge. In the case of several equally important sources of deformation along different axes, extra terms must be added to (24) and the analysis would become more complex.

To proceed it is necessary to use equation (24) to form the angular momentum $`𝑱`$ of the body. However as we are not modelling a rigid body, we must take care to allow for relative motion of one part with respect to another. Following  we will write the velocity of some point in the body as the sum of a rotational velocity with angular velocity $`𝛀`$ and a small velocity $`𝐮`$ relative to this rotating frame. We will call the frame that rotates at $`𝛀`$ the *body frame* although it is only in the rigid body limit that the body’s shape is fixed with respect to this frame. In other words the velocity of some particle making up the body is the sum of the body frame velocity $`𝛀\times 𝐫`$ at that point $`𝐫`$ plus the velocity $`𝐮`$ of the point relative to the body frame. Then

$$J_a=I_{ab}\mathrm{\Omega }_b+h_a,$$
(25)
where the possibly time-varying moment of inertia is defined in the usual way:

$$I_{ab}=_V\rho (x_cx_c\delta _{ab}x_ax_b)𝑑V,$$
(26)
while $`h_a`$ is the angular momentum of the body *relative* to this frame:

$$h_a=_V\rho ϵ_{abc}x_bu_c𝑑V.$$
(27)
We will neglect the $`h_i`$ term when constructing a free precessional motion, as it can be shown that $`h_i`$ is small in a well-defined sense . Therefore we will simply write

$$J_a=I_{ab}\mathrm{\Omega }_b$$
(28)
Having formulated the problem in this manner it is straightforward to show that the free precession of an elastic body is similar to that of a rigid one. First write down the angular momentum using (24) and (28). Referring all of our tensors to the body frame, with the 3-axis along $`𝒏_𝒅`$:

$$𝑱=(I_{0,S}+2/3\mathrm{\Delta }I_\mathrm{\Omega }1/3\mathrm{\Delta }I_d)𝛀+\mathrm{\Delta }I_d\mathrm{\Omega }{}_{3}{}^{}𝒏_{𝒅}^{}.$$
(29)
This shows that $`𝑱`$, $`𝛀`$ and $`𝒏_𝒅`$ are coplanar. As the angular momentum is constant this plane must rotate about $`𝑱`$. As in the rigid body case, we will refer to this as the *reference plane*. See figure 2. Taking the components of (29) we obtain:

$$J_1=(I_{0,S}+2/3\mathrm{\Delta }I_\mathrm{\Omega }1/3\mathrm{\Delta }I_d)\mathrm{\Omega }_1I_1\mathrm{\Omega }_1,$$
(30)
$$J_2=(I_{0,S}+2/3\mathrm{\Delta }I_\mathrm{\Omega }1/3\mathrm{\Delta }I_d)\mathrm{\Omega }_2I_1\mathrm{\Omega }_2,$$
(31)
$$J_3=(I_{0,S}+2/3\mathrm{\Delta }I_\mathrm{\Omega }+2/3\mathrm{\Delta }I_d)\mathrm{\Omega }_3I_3\mathrm{\Omega }_3.$$
(32)
These equations show that despite the triaxiality of $`𝑰`$ the angular momentum components themselves are structurally equivalent to those of a rigid symmetric top. The equations of motion of the body (i.e. Euler’s equations) involve only the components of $`𝑱`$ and $`𝛀`$. Therefore equations (30)—(32) show that the free precession of the triaxial body is formally equivalent to that of a rigid symmetric top. We can think of the elastic body as having an *effective* moment of inertia tensor $`diag[I_1,I_1,I_3]`$. Note that the *effective oblateness* $`I_3I_1`$ is equal to $`\mathrm{\Delta }I_d`$.

Now introduce standard Euler angles to describe the body’s orientation, with the polar axis along $`𝑱`$. Let $`\theta `$ and $`\varphi `$ denote the polar and azimuthal coordinates of the deformation axis, while $`\psi `$ represents a rotation about this axis. We refer to $`\theta `$ as the *wobble angle*. Taking the ratio of components $`J_1`$ and $`J_3`$ using (30) and (32) at an instant when $`\mathrm{\Omega }_2=0`$ we obtain

$$\mathrm{tan}\gamma =\frac{I_3}{I_1}\mathrm{tan}\theta ,$$
(33)
where $`\gamma `$ denotes the $`(𝛀,𝒏_𝒅)`$ angle. See figure 2.

We will label the angle between $`𝑱`$ and $`𝛀`$ as $`\widehat{\theta }`$:

$$\widehat{\theta }=\gamma \theta .$$
(34)
This angle is much smaller than $`\theta `$, as can be seen by linearising (33) in $`\mathrm{\Delta }I_\mathrm{\Omega }`$ and $`\mathrm{\Delta }I_d`$ to give

$$\widehat{\theta }=\frac{\mathrm{\Delta }I_d}{I_3}\mathrm{sin}\theta \mathrm{cos}\theta .$$
(35)
Note that according to our conventions, when the deformation bulge is oblate $`\mathrm{\Delta }I_d`$ and $`\widehat{\theta }`$ are positive, but when the deformation bulge is prolate $`\mathrm{\Delta }I_d`$ and $`\widehat{\theta }`$ are negative.

We can decompose the angular velocity according to

$$𝛀=\dot{\varphi }𝒏_𝑱+\dot{\psi }𝒏_𝒅.$$
(36)
Substituting this into equation (29) and resolving along $`𝒏_𝑱`$ and $`𝒏_𝒅`$ gives

$$J=I_1\dot{\varphi },$$
(37)
$$\dot{\psi }=\frac{\mathrm{\Delta }I_d}{I_1}\mathrm{\Omega }_3,$$
(38)
where $`J`$ denotes the magnitude of the angular momentum. Note that when $`\mathrm{\Delta }I_\mathrm{\Omega }=0`$ the above formulae reduce to the familiar rigid body equations.

Thus the motion is simple. As viewed from the inertial frame the deformation axis rotates at a rate $`\dot{\varphi }`$ in a cone of half-angle $`\theta `$ about the angular momentum vector. This angular velocity is sometimes called the *inertial precession frequency*. The centrifugal bulge rotates around the angular momentum vector also, but—for oblate deformations—on the opposite side of $`𝑱`$, making an angle $`\widehat{\theta }\gamma \theta `$ with $`𝑱`$. Superimposed upon this is a rotation about the deformation axis at a rate $`\dot{\psi }`$, known as the *body frame precession frequency* or sometimes simply the *precession frequency*. This frequency is negative for an oblate distortion and positive for a prolate one.

## V Radiation reaction for an elastic body: Energy and angular momentum balance

Here we derive the wobble damping time $`\tau _\theta `$ for elastic bodies, based on energy and angular momentum balance. Once fully underway, the derivation is just a couple lines. But to understand it, it is useful to carry along a simple, physical model for the deformed crust. (However our derivation will actually be completely general). Here is the model: take some non-rotating, spherical NS, and stretch a rubber band around some great circle on the crust. We shall refer to this great circle as the NS’s equator. Obviously the effect of the rubber band is to make the NS slightly prolate (but still axisymmetric). To get an oblate shape, you can instead imagine sewing compressed springs into the surface of the crust at the equator. For definiteness, let the potential energy of the band (or springs) be $`V=\frac{1}{2}ϵl^2`$, where $`l`$ is its length. So $`ϵ`$ is positive for the rubber band (prolate deformation, $`\mathrm{\Delta }I_d<0`$) and negative for the springs (oblate deformation, $`\mathrm{\Delta }I_d>0`$). Now give the NS angular momentum $`𝑱`$ about some axis that is not quite perpendicular to the equator. We now have our deformed, wobbling NS. We consider the equation of state of the star and the value $`ϵ`$ to be fixed once and for all, and consider how the energy of the system (star + band) varies as a function of its total angular momentum $`J`$ and the wobble angle $`\theta `$ (the angle between $`𝑱`$ and the perpendicular to the equator); i.e., we consider $`E(J,\theta )`$. We will be concerned with small wobble angle, so let us expand $`E(J,\theta )`$ as a Taylor series in $`J`$ and $`\theta `$:

$$E(J,\theta )=E_0+\frac{1}{2}BJ^2+\frac{1}{24}CJ^4+\frac{1}{2}Fϵ\theta ^2J^2+\mathrm{}$$
(39)
Here $`E_0`$ is defined to be the energy of the (star + band) at zero $`J`$, and $`B,C`$, and $`F`$ are some expansion coefficients that in principle depend on the physical properties of the (star + band). Fortunately we will soon see that there are simple relations between $`B,C`$, and $`F`$ and previously-defined physical parameters, such as $`\mathrm{\Delta }I_d`$. Our ultimate goal is to obtain the two partial derivatives on the right-hand side of equation (5), where $`E`$ now denotes the *total* energy.

First, to see that no lower order terms (such as $`J,\theta J,\theta ^2`$, or $`\theta J^2`$ terms) can appear in the expansion (39), note that the $`J=0`$ configuration corresponds to the minimum of the potential energy of the (star + band) system. Displacements of the (star + band) are first order in $`J^2`$, so changes in the potential energy of (star + band) are $`𝒪(J^4)`$. Thus terms in $`E(J,\theta )`$ that are $`J^2`$ are kinetic energy pieces. These terms with a $`J^2`$ in them are clearly just $`\frac{1}{2}(I_0^1)^{ab}J_aJ_b`$, where $`I_0^{ab}`$ is defined to be the inertia tensor of the (star + band) at $`J=0`$. (Corrections to the star’s $`I^{ab}`$ first enter the energy at $`𝒪(J^4)`$.) We write $`I_0^{ab}`$ as

$$I_0^{ab}=I_{0,S}\delta ^{ab}+\mathrm{\Delta }I_d(n_d^an_d^b\frac{1}{3}\delta ^{ab}),$$
(40)
where $`I_{0,S}`$ represents the ‘spherical part’ of $`I_0^{ab}`$. Then

$$(I_0^1)^{ab}=\frac{1}{I_{0,S}}[\delta ^{ab}\left(\frac{\mathrm{\Delta }I_d}{I_{0,S}}\right)(n_d^an_d^b\frac{1}{3}\delta ^{ab}\left)\right]$$
(41)
where a term of $`𝒪(\mathrm{\Delta }I_d^2)`$ has been neglected. The kinetic energy part of $`E`$ is \[up to terms of $`𝒪(\mathrm{\Delta }I_d^2)`$ and $`𝒪(J^4)`$ \]

$$E_{kin}=\frac{J^2}{2I_{0,S}}\left[1\left(\frac{\mathrm{\Delta }I_d}{I_{0,S}}\right)\left(\frac{2}{3}\theta ^2\right)\right],$$
(42)
where we have used the small wobble angle result $`J_an_d^a=J(1\frac{1}{2}\theta ^2)`$. From Eq. (42) we immediately read off the values of $`B`$ and $`Fϵ`$ in expansion (39):

$`B=I_{0,S}^1\left[1{\displaystyle \frac{2}{3}}{\displaystyle \frac{\mathrm{\Delta }I_d}{I_{0,S}}}\right],`$ (43)
$`Fϵ=\mathrm{\Delta }I_d/(I_{0,S})^2,`$ (44)
and obtain the partial derivative

$$\frac{E}{\theta }|_J=J^2\theta \frac{\mathrm{\Delta }I_d}{(I_{0,S})^2}.$$
(45)
To compute the partial derivative in the numerator of equation (5) it is sufficient to consider the $`\theta 0`$ limit  so that

$$\mathrm{\Omega }=\frac{dE}{dJ}|_{\theta =0}=BJ+\frac{1}{6}CJ^3,$$
(46)
where $`\mathrm{\Omega }`$ denotes the spin frequency in the axisymmetric limit. It is related to the inertial precession frequency by

$$\mathrm{\Omega }=\dot{\varphi }(1\mathrm{\Delta }I_d/I_{0,S}).$$
(47)
The final physics inputs we need are

$$\dot{E}=\frac{2G}{5c^5}\frac{\mathrm{\Delta }I_d^2}{I_{0,S}}\dot{\varphi }^6\theta ^2$$
(48)
$$\dot{E}=\dot{\varphi }\dot{J}.$$
(49)
Eqs. (48) and (49) follows from the quadrupole formalism in the same way as for the rigid body.

The necessary pieces have been gathered; substituting into equation (5) gives

$`\dot{\theta }`$ $`=`$ $`{\displaystyle \frac{\dot{J}I_{0,S}^2}{\theta J^2}}{\displaystyle \frac{(\dot{\varphi }\mathrm{\Omega })}{\mathrm{\Delta }I_d}}`$ (50)
$`=`$ $`{\displaystyle \frac{2G}{5c^5}}({\displaystyle \frac{\mathrm{\Delta }I_d}{I_{0,S}}})^2I_{0,S}\dot{\varphi }^4\theta .`$ (51)
This is simply the same spin-down rate as for a rigid body, with the replacement $`(\mathrm{\Delta }I/I_1)ϵ_d`$. This is much longer than the timescale claimed by Bertotti and Anile  by a factor $`\mathrm{\Delta }I_\mathrm{\Omega }/\mathrm{\Delta }I_d`$, which is typically $`10^5`$ or higher.

Finally, the spin-down rate $`\ddot{\varphi }`$ can be obtained in the same way as for a rigid body, i.e. by differentiating $`\dot{\varphi }=J/I_1`$ and using equations (48) and (49). Strictly there will also be a term in $`\dot{I_1}`$, but this correction will be down by a factor of order $`(\mathrm{\Omega }/\mathrm{\Omega }_{max})^2`$. We then obtain the same spin-down as for a rigid body, again with the replacement $`\mathrm{\Delta }I\mathrm{\Delta }I_d`$:

$$\ddot{\varphi }=\frac{2G}{5c^5}\frac{\mathrm{\Delta }I_d^2}{I_{0,S}}\dot{\varphi }^5\theta ^2.$$
(52)
## VI Radiation reaction for an elastic body: Local force

We now give a second derivation of the wobble damping rate for an elastic star, by directly adding the gw radiation reaction force to the Newtonian equations of motion. Besides being a satisfying consistency check on the calculation in §4, by doing this second derivation correctly we can show where Bertotti and Anile  went astray.

As was the case for the rigid body, the Burke-Thorne potential will exert a torque on the spinning star. However, this is not the only effect of the radiation reaction force: It will distort the shape of the NS and thus its moment of inertia. The equation describing the precession is then of the from

$$\frac{d}{dt}[(𝑰_𝑵+\delta 𝑰_{𝑩𝑻})𝛀]=𝑻,$$
(53)
where $`𝑰_𝑵`$ denotes the Newtonian part of the moment of inertia tensor, $`\delta 𝑰_{𝑩𝑻}`$ the perturbation in this tensor due to the Burke-Thorne force, and $`𝑻`$ the Burke-Thorne torque. It was the $`\delta 𝑰_{𝑩𝑻}`$ terms that were not included by Bertotti and Anile. Fortunately, these can also be calculated explicitly, as we show below.

### A Effect of $`\mathrm{\Phi }^{RR}`$ on the NS’s Shape

It is perhaps surprising that one can explicitly determine the effect of $`\mathrm{\Phi }^{RR}`$ on the NS’s moment of inertia, since the answer would seem to depend on the NS’s mass and the details of its equation of state; i.e., one might worry that extra parameters must be specified even to make the problem well-defined. However the point is that (from symmetry arguments) the perturbation $`\mathrm{\Delta }I_{ij}`$ depends only on a single physical parameter, and this parameter already appears in our Newtonian equations of motion. That parameter is $`\mathrm{\Delta }I_\mathrm{\Omega }/\mathrm{\Omega }^2`$, the amount of oblateness caused “per unit centrifugal force”.

The point is that both the centrifugal and radiation reaction forces have the very special property that they grow linearly with distance from the center of the star. This fact, coupled with symmetry arguments, is enough to determine $`\mathrm{\Delta }I_{ij}`$ in terms of $`\mathrm{\Delta }I_\mathrm{\Omega }/\mathrm{\Omega }^2`$; no new physical parameters have to be introduced.

Let $`\mathrm{\Phi }^\mathrm{\Lambda }`$ be some external potential of the form $`\mathrm{\Phi }^\mathrm{\Lambda }\mathrm{\Lambda }^{ab}x_ax_b`$, where $`\mathrm{\Lambda }_{ab}`$ is some trace-free tensor. Allow this potential to act on the non-rotating (and so spherically symmetric) NS; it will induce a perturbation $`\mathrm{\Delta }I^{ab}`$ in the NS’s inertia tensor. Since the background is spherically symmetric, the only possibility (to first order in the perturbation) is that $`\mathrm{\Delta }I^{ab}=C\mathrm{\Lambda }^{ab}`$, where $`C`$ is some constant (i.e, independent of $`\mathrm{\Lambda }^{ab}`$).

We can determine $`C`$ as follows. Decompose the centrifugal potential into a spherically symmetric and a trace-free piece:

$$\frac{1}{2}\mathrm{\Omega }^2(\delta ^{ab}n_\mathrm{\Omega }^an_\mathrm{\Omega }^b)x_ax_b=\frac{1}{3}\mathrm{\Omega }^2x^2+\mathrm{\Lambda }_\mathrm{\Omega }^{ab}x_ax_b,$$
(54)
where $`\mathrm{\Lambda }_\mathrm{\Omega }^{ab}=\frac{1}{2}\mathrm{\Omega }^2(n_\mathrm{\Omega }^an_\mathrm{\Omega }^b\frac{1}{3}\delta ^{ab})`$. For small $`\mathrm{\Omega }`$ the perturbed inertia tensor is $`\mathrm{\Delta }I^{ab}=\mathrm{\Delta }I_\mathrm{\Omega }(n_\mathrm{\Omega }^an_\mathrm{\Omega }^b\frac{1}{3}\delta ^{ab})`$, so the constant $`C`$ is just $`2\mathrm{\Delta }I_\mathrm{\Omega }/\mathrm{\Omega }^2`$.

The radiation reaction potential for the freely precessing elastic body can be found by substituting the radiation reaction free motion into equation (14) to give:

$$\mathrm{\Phi }^{RR}=\frac{G}{5c^5}x^ax^b\left[\mathrm{\Delta }I_d\frac{d^5}{dt^5}(𝒏_{𝒅}^{}{}_{a}{}^{}𝒏_{𝒅}^{}{}_{b}{}^{})+\mathrm{\Delta }I_\mathrm{\Omega }\frac{d^5}{dt^5}(𝒏_{𝛀}^{}{}_{a}{}^{}𝒏_{𝛀}^{}{}_{b}{}^{})\right].$$
(55)
The first term is the potential caused by the motion of the deformation bulge, the second by the centrifugal bulge. The differentiations of the unit vectors are straightforward. In the case where $`\theta 1`$ we can approximate $`𝒏_𝒅𝒏_𝑱+\theta 𝒏_\mathbf{}𝑱`$ and $`𝒏_𝛀𝒏_𝑱\widehat{\theta }𝒏_\mathbf{}𝑱`$, where $`𝒏_\mathbf{}𝑱`$ is the unit vector in the reference plane which lies perpendicular to $`𝐉`$ and points towards $`𝒏_𝒅`$. We then find

$$\mathrm{\Phi }^{RR}=\frac{G}{5c^5}\dot{\varphi }^5x^ax^b\dot{\varphi }^5[\mathrm{\Delta }I_d\theta \mathrm{\Delta }I_\mathrm{\Omega }\widehat{\theta }](\widehat{𝒗}_a𝒏_{𝑱}^{}{}_{b}{}^{}+𝒏_{𝑱}^{}{}_{a}{}^{}\widehat{𝒗}_b).$$
(56)
Here $`\widehat{𝒗}`$ denotes a unit vector $`𝒏_𝑱\times 𝒏_\mathbf{}𝑱`$. Using the prescription described above, these radiation reaction potentials can be converted immediately into perturbations of the moment of inertia tensor:

$$\delta 𝑰_{𝑩𝑻}=\frac{2G}{5c^5}\dot{\varphi }^3[\mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }\theta (\mathrm{\Delta }I_\mathrm{\Omega })^2\widehat{\theta }](\widehat{𝒗}𝒏_𝑱+𝒏_𝑱\widehat{𝒗}).$$
(57)
### B Adding $`\mathrm{\Phi }^{RR}`$ to Equations of Motion

It now remains to compute the torque $`𝑻`$ using equation (18). We obtain four terms, corresponding to the expansion of the product of $`𝐈\text{ }`$ with its fifth time derivative. Again linearising with respect to $`\theta `$ we obtain

$$𝑻=\frac{2G}{5c^5}\dot{\varphi }^5[\mathrm{\Delta }I_d^2\theta \mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }\widehat{\theta }+\mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }\theta \mathrm{\Delta }I_\mathrm{\Omega }^2\widehat{\theta }]𝒏_\mathbf{}𝑱.$$
(58)
Define $`ϵ_\mathrm{\Omega }\mathrm{\Delta }I_\mathrm{\Omega }/I_{0,S}`$ and $`ϵ_d\mathrm{\Delta }I_d/I_{0,S}`$ Note our definition of $`ϵ_\mathrm{\Omega }`$ differs by a factor $`2/3`$ from , who set $`ϵ_\mathrm{\Omega }\frac{2}{3}\mathrm{\Delta }I_\mathrm{\Omega }/I_{0,S}`$. Then the terms on the rhs of (58) stand in the ratio $`ϵ_d/ϵ_\mathrm{\Omega }:ϵ_d:1:ϵ_\mathrm{\Omega }`$. We are now in a position to write down the equation for $`d(𝑰_N𝛀)/dt`$. Using Eq. (57) and the Newtonian motion to compute $`d[(\delta 𝑰_{𝑩𝑻})𝛀]/dt`$, and neglecting terms of order $`\theta ^2`$, we find that Eq. (53) reduces to

$$\frac{d}{dt}(𝑰_N𝛀)+\frac{2G}{5c^5}\dot{\varphi }^5\theta [\mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }\theta (\mathrm{\Delta }I_\mathrm{\Omega })^2\widehat{\theta }]𝒏_\mathbf{}𝑱$$
(59)
$`={\displaystyle \frac{2G}{5c^5}}\dot{\varphi }^5[\mathrm{\Delta }I_d^2\theta \mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }\widehat{\theta }+\mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }\theta \mathrm{\Delta }I_\mathrm{\Omega }^2\widehat{\theta }]𝒏_\mathbf{}𝑱.`$We see that the last two terms on the rhs are cancelled by terms on the lhs. This leaves

$$\frac{d}{dt}(𝑰_N𝛀)=\frac{2G}{5c^5}\dot{\varphi }^5[\mathrm{\Delta }I_d^2\theta \mathrm{\Delta }I_d\mathrm{\Delta }I_\mathrm{\Omega }\widehat{\theta }]𝒏_\mathbf{}𝑱.$$
(60)
The problem has reduced to a rigid-body Newtonian one, with the two torque terms indicated on the right-hand side. The terms stand in the ratio $`1:ϵ_\mathrm{\Omega }`$. In fact, the dominant term is the same as that obtained in the rigid body case with the change $`\mathrm{\Delta }I\mathrm{\Delta }I_d`$.

We therefore find that the alignment rate as calculated using the local Burke-Thorne formalism agrees with the flux-at-infinity method. The previous force-based calculation of Bertotti and Anile  failed to include the deformation $`\delta 𝑰_{𝑩𝑻}`$, so that the cancellations in equation (59) described above did not occur.

Finally, it is easy to show that even when the approximations $`\theta 1`$, $`ϵ_d1`$ are *not* employed, the effective torques due to the $`\delta 𝑰_{𝑩𝑻}`$ terms are still perpendicular to $`𝐉`$, so the spin-down $`\ddot{\varphi }`$ using this local formalism is necessarily the same as in the flux-at-infinity method.

## VII Allowance for a liquid core

We have successfully described the effects of gravitational radiation reaction on an elastic precessing body. We will now briefly describe how to extend this result to the realistic case where the star consists of an elastic shell (the crust) containing a liquid core. The Earth itself is just such a body, and the form of its free precession was considered long ago. We will base our treatment on that of Lamb , who considered a rigid shell containing an incompressible liquid of uniform density. To make the problem tractable the motion of the fluid was taken to be one of uniform vorticity. We will assume the ellipticity of the shell, and also the ellipticity of the cavity in which the fluid resides, are small. Then the small angle free precession of the combined system can be found by means of a normal mode analysis of the equations of motion .

The key points are as follows: The fluid’s angular velocity vector does not significantly participate in the free precession. Instead it remains pointing along the system’s total angular momentum vector. The shell precesses about this axis in a cone of constant half-angle. The fluid exerts a force on the shell such that the shell’s body frame precession frequency is increased in magnitude, so that:

$$\dot{\psi }=\dot{\varphi }\frac{\mathrm{\Delta }I}{I_{\mathrm{crust}}}$$
(61)
where $`\mathrm{\Delta }I`$ denotes the difference between the 1 and 3 principal moments of inertia of the *whole body*, not just the shell.

We now wish to calculate the alignment rate of such a body due to gravitational radiation reaction. The averaged energy and angular momentum fluxes, as well as the instantaneous torque, depend only upon the orientation of the mass quadrupole of the body, and so are exactly the same as if the body were rigid, i.e. equations (1), (2) and (19) apply. Equations giving the kinetic energy and angular momentum of the body are given in Lamb . These can be used to obtain the partial derivatives that appear in equation (5). Explicitly, we find

$$\frac{E}{J}|_\theta =\mathrm{\Omega }=\dot{\varphi }+\dot{\psi }$$
(62)
and

$$\frac{E}{\theta }|_J=\dot{\varphi }^2\theta \mathrm{\Delta }I.$$
(63)
(See Jones  for a detailed derivation.)

These lead to an alignment timescale that is $`I_{\mathrm{crust}}/I`$ shorter than that of equation (13). This result is confirmed using the local torque formulation, where

$$\dot{\theta }=\frac{T_J}{I_{\mathrm{crust}}}.$$
(64)
In the realistic case where both crustal elasticity and core fluidity are taken into account we can combine the above arguments as described by Smith and Dahlen , i.e. we can take the rigid result and put $`II_{\mathrm{crust}}`$ and $`\mathrm{\Delta }I\mathrm{\Delta }I_d`$ to give

$`\dot{\psi }=\dot{\varphi }{\displaystyle \frac{\mathrm{\Delta }I_d}{I_{\mathrm{crust}}}}`$ (65)
$`\dot{\theta }={\displaystyle \frac{2G}{5c^5}}{\displaystyle \frac{\mathrm{\Delta }I_d^2}{I_{\mathrm{crust}}}}\dot{\varphi }^4\theta .`$ (66)
## VIII Conclusions

We have shown that the gw damping time for wobble in realistic NS’s has the same form as for rigid bodies, but with the replacement $`\mathrm{\Delta }I^2/I_1\mathrm{\Delta }I_d^2/I_{\mathrm{crust}}`$. This given an alignment timescale of:

$$\tau _\theta =1.8\times 10^5\mathrm{yr}\left(\frac{I_{\mathrm{crust}}}{10^{44}\mathrm{g}\mathrm{cm}^2}\right)\left(\frac{10^{38}\mathrm{g}\mathrm{cm}^2}{\mathrm{\Delta }I_d}\right)\left(\frac{\mathrm{kHz}}{\nu _s}\right)^4.$$
(67)
For the Crab, taking $`ϵ_d3\times 10^9`$, this gives $`\tau _\theta 5\times 10^{13}`$yrs–much longer than the age of the universe. For an accreting NS with $`ϵ_d10^7`$ and $`\nu _s300`$ Hz, we estimate $`\tau _\theta 2\times 10^8`$ yrs.

Our basic conclusion, then, is that gw backreaction is sufficiently weak that other sources of dissipation probably dominate. Unfortunately, even for the Earth the dissipation mechanisms are not well understood . Early estimates of Chau and Henriksen , which considered dissipation within the neutron star crust, suggested that wobble would be damped in around $`10^6`$ free precession periods, i.e. over a time interval of $`10^6/(ϵ_d\nu _s)`$. A more recent study of Alpar and Sauls  argued that the dominant dissipation mechanism will be due to imperfect coupling between the crust and the superfluid core. They estimate that the free precession will be damped in (at most) $`10^4`$ free precession periods. In contrast, according to equation (67), the gw damping time is in excess of $`10^8\left(\frac{\mathrm{kHz}}{\nu _s}\right)^3`$ free precession periods. On the basis of these estimates, it seems likely that internal damping will dominate over gravitational radiation reaction in *all* neutron stars of interest. Note however, that while internal dissipation damps wobble for *oblate* deformations, we expect that internal dissipation causes the wobble angle to *increase* in the prolate ($`\mathrm{\Delta }I_d<0`$) case.

A study of the gravitational wave detectability of realistic neutron stars undergoing free precession, including a discussion of other astrophysical mechanisms which might affect the evolution of the motion, will be presented elsewhere (Jones, Schutz and Andersson, in preparation).

###### Acknowledgements.

We thank N. Andersson and B. F. Schutz for discussions. This research was supported by NASA grant NAG5-4093 and PPARC grant PPA/G/1998/00606.


text:
# Near-infrared and Br𝛾 observations of post-AGB starsBased on observations made at the European Southern Observatory, La Silla, Chile; the Australia Telescope Compact Array, which is funded by the Commonwealth of Australia for operations as a National Facility managed by CSIRO; Siding Spring Observatory.

## 1 Introduction

A most intriguing challenge is to understand how Asymptotic Giant Branch (AGB) stars transform their surrounding mass-loss shells in a couple of thousand years into the variety of shapes and sizes observed in Planetary Nebulae (PNe). There are a number of theories currently being investigated. In the generalized interacting stellar wind model, a variety of axisymmetric PN shapes are obtained by the interaction of a very fast central star wind with the progenitor AGB circumstellar envelope (Kwok Kwok82 (1982)), when the latter is denser near the equator than the poles (Frank et al. Frank93 (1993) ). Sahai & Trauger (Sahai98 (1998)) proposed that the primary agent for shaping PNe is not the density contrast, but a high speed collimated outflow of a few 100 km/s. No consensus about the dominant physical process responsible for the shaping of PNe has emerged so far, but there is agreement that they occur during the early AGB-to-PN transition stage. However, details of the rate of evolution, the strength of the stellar wind, and the impact of ionization across the transition phase are very poorly known. In order to test model predictions it is essential to study nebulae from the early post-AGB phase through the very early PN phases. However, post-AGB objects are difficult to find, because this phase is very short and the star is usually obscured by a thick circumstellar dust shell.

One way of identifying new post-AGB stars is through their dust emission. In our search for new obscured PNe and post-AGB stars, we selected candidates from the IRAS Point Source Catalogue based on infrared colors typical of PNe. When we observed these candidates in the radio continuum at 6 cm, on average 20 % of the objects were detected (Van de Steene & Pottasch VdSteene93 (1993), VdSteene95 (1995)). Subsequent optical spectroscopy showed that the PN candidates detected in the radio have emission line spectra typical of PNe (Van de Steene et al. 1996a , 1996b ). However, the question remained: what is the evolutionary status of those IRAS sources with colors typical of PNe, which had no detectable ionization in the radio (i.e. fluxes below 3 mJy)? It is possible that a few large, low surface-brightness PNe have escaped detection. Others of the remaining non-identified PN candidates could be very young and small PNe. Their radio flux would have been below our detection threshold and they wouldn’t be identified in H$`\alpha `$ due to extinction. However, the true evolutionary status of most objects with IRAS colors typical of PNe has remained unknown.

We have calculated evolutionary tracks for post-AGB stars in the IRAS color-color diagram (van Hoof et al. vHoof97 (1997)). From this study, it became clear that post-AGB objects and PNe could be located in the same region in the IRAS color-color diagram. No other types of objects seem to have typically these particular IRAS colors. We therefore adopted the working hypothesis that the non-detected PN candidates are AGB-to-PN transition objects. The goal of this project is to determine the evolutionary status of these post-AGB candidates, their physical properties, and to relate them to their morphology.

To ensure the identification of the correct counterpart of the IRAS source and obtain accurate positions for follow-up observations, we first imaged the post-AGB candidates in the N-band. In order to investigate whether the sources have some ionization, we obtained Br$`\gamma `$ spectra. We reobserved the sources showing Br$`\gamma `$ emission at 3 cm and 6 cm. To confirm the identification, obtain finding charts, and photometry, we took high resolution near-infrared (JHKL) images.

The next section describes the observations. In Section 3 we describe the data reduction and analysis, in Section 4 we discuss the individual objects and in Section 5 we discuss the general results. We give a summary in the last section.

## 2 Observations

### 2.1 Imaging

#### 2.1.1 N-band

Mid-infrared images in the N-band were obtained with the Thermal Infrared Multi Mode Instrument (timmi) at the European Southern Observatory 10–12 May 1995. The camera featured a 64 $`\times `$ 64 element Gallium doped Silicon array with good cosmetic quality and a quantum efficiency of 25 %. We used the N-band filter: $`\lambda _\mathrm{c}`$ = 9.70 $`\mu `$m and $`\mathrm{\Delta }\lambda `$ = 5.10 $`\mu `$m with a setup yielding a pixel size of 0$`\stackrel{}{.}`$66, and a FOV of $``$ 40″.

For the cancelation of the strong background radiation, timmi is operated in chopping and nodding mode. The chopping frequency is several hertz and nodding is done once or twice per minute. In chopping mode a pair of exposures is obtained: one exposure contains object plus sky and the second exposure, 20″ away, only sky. To obtain the images in nodding mode the telescope was moved by an angle exactly matching the chopper throw on the sky and the observation in chopping mode was repeated. Because this observation is exactly 180 degrees out of phase with the first one, the observed object shows up as a negative measurement. Subtracting both images resulted in a clean image corrected for first and second order effects introduced by the strong and inhomogeneous background radiation. To eliminate erratic pixels, the multiple images of each object were median filtered with a high threshold. Finally the images were Fourier transformed and the midas standard low bandpass Butterworth filter was applied to better reveal the sources.

A log of the observations can be found in Table 1. The weather was poor during both nights. Consequently, only the brighter IRAS sources could be detected, while under favorable weather conditions timmi would be ten times more sensitive than IRAS. All but one of the undetected sources have an IRAS 12-$`\mu `$m flux below 4 Jy. Six sources having a 12-$`\mu `$m flux below 4 Jy were detected when the cloud conditions improved temporarily. The other 14 detected sources have a larger IRAS 12-$`\mu `$m flux. We think that non-detections are due to the cloudy weather rather than faulty IRAS positions or inaccurate telescope pointing.

#### 2.1.2 Near-infrared images

JHKL-band images of the 20 post-AGB candidates detected in the N-band were obtained with caspir at the 2.3-m telescope at Siding Spring Observatory in Australia on 11 and 12 March 1998. The detector is a SBRC CRC463 $`256\times 256`$ InSb array which is sensitive from approximately 0.9 $`\mu `$m to 5.5 $`\mu `$m. The pixel size was 0$`\stackrel{}{.}`$25 and the total field of view 1′. Each object was observed 60 times for 0.3 s at 12″ north and 12″ south of its nominal position. The objects were observed in order of right ascension: IRAS 10256$``$5628 to 15144$``$5812 on 11 March 1998 and IRAS 15544$``$5332 to 17088$``$4221 on 12 March 1998. Both nights were photometric. Standard stars were observed at the beginning of the observations each night. The seeing was sub-arcsecond.

### 2.2 Observing the Br$`\gamma `$ line

Using the improved positions from the timmi images, we obtained infrared spectra with irspec at the NTT at ESO 13 and 14 May 1995.

The spectra are centered at Br$`\gamma `$ and have a resolution of $`\lambda /\mathrm{\Delta }\lambda 2000`$. The SBRC $`62\times 58`$ array gives a total wavelength coverage from 2.15 $`\mu `$m to 2.18 $`\mu `$m. The resolution in the spatial direction is 2$`\stackrel{}{.}`$26/pixel. Since the outer edges of the array are vignetted, the total field of view is approximately 1$`\stackrel{}{.}`$5. The slit width is 2 pixels or 4$`\stackrel{}{.}`$51. The slit orientation was north-south. First we pointed the telescope to a bright star, and subsequently offset the telescope to the timmi position. The pointing accuracy of the NTT has an rms of about 1″. Subsequently we moved the slit east- and westwards a few arcseconds, to ‘peak-up’ the strongest signal in order to center the object in the slit. The resulting offsets were never larger than 6″.

To subtract the sky emission the beam-switching technique was used. Because all but one (IRAS 13428$``$6232) of the sources are much smaller than 20″, two spectra were taken with the source at different positions on the array, such that when the exposures were subtracted, the spectra did not overlap. For IRAS 13428$``$6232 one integration was taken on source and one on the sky and these images were subsequently subtracted

The weather was very good, and the humidity was around 15 %. Within the accuracy limits of spectro-photometry the night was photometric, even if the conditions slightly deteriorated at the end of the night. Every 30 to 60 min a standard star was observed. A complete log of the observations can be found in Table 2. Note that the observations include four PNe from Van de Steene & Pottasch (VdSteene95 (1995)).

### 2.3 Radio continuum observations

The 6 objects showing Br$`\gamma `$ in emission were observed with the Australia Telescope Compact Array for 12 h on 14 February 1997 and for 2 h on 15 February 1997. The array configuration was $`\mathrm{\#}`$6A with the antennae at stations 3, 11, 16, 30, 34, and 37. The shortest baseline was 337 m and the longest 5939 m. The bandwidth was 128 MHz divided in 32 channels centered at 4800 MHz and 8640 MHz corresponding to 6 cm and 3 cm respectively. The sources were observed simultaneously at both frequencies. To obtain adequate spatial coverage we cycled through the sample of 6 IRAS objects and 3 calibrators once an hour: 1338$``$59 C, 13428$``$6232, 14488$``$5405, 1511$``$55 C, 15144$``$5812, 15544$``$5332, 1600$``$48 C, 16594$``$4656, 17009$``$4154 (with C denoting the calibrators). Each object was observed for 8 min and every phase calibrator for 3 min. To avoid artifacts in the field center we used an offset of 30″ in declination. At the beginning and the end of each shift the primary flux density calibrator 1934$``$638 was observed for 5 min to 10 min. The weather was unstable and humid.

## 3 Data reduction and analysis

### 3.1 Imaging

#### 3.1.1 N-band

At 10 $`\mu `$m, fields are not crowded: for each of our images we had only one object within a relatively small field of view of 40″. This leaves no doubt about the identification of the correct counterpart in the timmi field at this wavelength, even in non-photometric weather conditions. The disadvantage is that there are no other points of reference to improve upon the IRAS coordinates other than the telescope pointing.

The positions and Full Width at Half Maximum (FWHM) of the sources were determined from the images after filtering. Because of the variable sky conditions, the background level around the sources was badly defined. Consequently the determination of the (standard) star’s PSF was very uncertain. To remedy this, we subtracted the background in the following way: we made an image in which we replaced a circular area surrounding the source by a flat surface fitted to the background outside this circle and subtracted this image from the original. The resulting frame contained only the object against a virtual zero background. Next a two dimensional Gaussian was fitted to the source to determine the position and FWHM more accurately.

At regular intervals during the two nights, we pointed the telescope to a standard star and put it exactly on the cross hairs. Next, we took an image to determine exactly the corresponding position on the array. Similarly, for each object, we calibrated the telescope pointing using the closest bright SAO star, before doing an offset to the IRAS position. There was always only one source in the field, which left no doubt about the identification of the IRAS source in the N-band image. We assumed that the IRAS position corresponds to the reference position on the array as determined from the standard stars. We measured the offset from this array position to the position of the source in the timmi image, and adopted this as the improved position of the IRAS source. The difference between the IRAS and timmi positions were mostly less than 10″ which is in agreement with what was found in positional difference between the radio detections and the IRAS positions (Van de Steene & Pottasch VdSteene93 (1993), VdSteene95 (1995)). The uncertainty in the timmi positions was estimated to be less than $``$5″ and is due to the fact that in these regions with high extinction, the nearest SAO star could be several arcminutes away and the pointing of the ESO 3.6-m was not very good.

The median seeing was 2$`\stackrel{}{.}`$7 FWHM during the first night and, even worse, 3$`\stackrel{}{.}`$4 FWHM, during the second. Fifteen objects remained unresolved. Five appear to be extended at 10 $`\mu `$m. The contour plots of the extended sources are presented in Fig. 8. The morphology of the resolved objects is similar to young PNe. The size of the extended sources is given in Table 3. One showed Br$`\gamma `$ in absorption, two in emission, and in two we didn’t detect any Br$`\gamma `$ in absorption or emission.

#### 3.1.2 Near-infrared

The JHKL-band images were reduced in iraf using standard procedures as described in the caspir manual (McGregor McGregor94 (1994)). The two images taken 24″ apart in declination were combined. Because the objects were most prominent in the K-band, the resulting K-band images are presented in Fig. 9.

The positions of the objects determined from these images are listed in Table 4. We used the imcoor package in iraf and the positions of reference stars from the United States Naval Observatory Catalogue (USNO release A1.0, available via the ESO skycat tool) to determine an accurate position. If the object had a counterpart in the catalogue, we adopted the USNO catalogue position. The positional uncertainty is similar to the uncertainty in the USNO positions, about 1″. The object for which we have obtained the Br$`\gamma `$ spectrum is indicated with a box in Fig. 9. Note that in some fields more than one star was in the slit. The IRAS position is indicated with a cross in Fig. 9.

Although the K-band seeing was sub-arcsecond, most of the objects are unresolved at this level. The sizes of the extended objects are presented in Table 3. All objects which are extended in the K-band show Br$`\gamma `$ in emission, but not vice versa. Obviously, comparison of morphology in the N\- and K-band is not a good identification tool.

The photometry was done with the program qphot in iraf. The JHKL-band magnitudes of the objects were determined for each of the two images separately, and their agreement checked before averaging. The average difference between the two measurements was 0.01 mag for both nights in JHK and 0.05 mag for the L-band. The near-infrared magnitudes are presented in Table 5. The uncertainties are estimated at 0.05 mag in JHK and 0.1 mag in the L-band, including measurements errors and the uncertainty in the correction for atmospheric extinction. Table 5 also contains estimates for the K magnitudes derived from the Br$`\gamma `$ spectra (see also Sect. 3.2). Usually these magnitudes are a bit fainter than the ones determined from the caspir images. This is probably due to slit loss.

Table 5 also presents near-infrared photometry from the literature, when available. The literature values are based on aperture photometry rather than imaging. In crowded regions they are often brighter, likely due to contamination by neighboring objects. In all cases where the literature value was very different from our magnitude, the two values were not associated with the same star, and we were able to identify the star measured by the other authors in our images.

Near-infrared imaging is preferred to aperture photometry to identify the object in the field. It can provide very accurate coordinates and finding charts for follow-up observations. Given photometric weather and good seeing conditions, the magnitudes won’t be contaminated by neighboring stars. When several near-infrared bands are available, colors can be used as a secondary means of identification. The counterpart often is the reddest, but not necessarily the brightest object in the near infrared. Due to their thick circumstellar dust shells, post-AGB stars often are very reddened.

### 3.2 The Br$`\gamma `$ spectra

The data reduction was done using the standard reduction macros contained in the midas image processing system. After flatfielding, sky subtraction, and rectifying the spectrum the resultant image contained a positive and a negative spectrum of the source. The next steps were to extract the positive and negative beam, invert the negative beam, and average the two. The spectra were accurately wavelength calibrated using the sky-emission lines present in the non-subtracted images. The rest wavelength of these lines were taken from Oliva & Origlia (Oliva92 (1992)). Using the standard star spectra, a flux conversion factor was determined to calibrate all spectra.

The values of the continuum flux at Br$`\gamma `$ (2.166 $`\mu `$m) have been converted to Johnson K-band magnitudes for all sources and are listed in Table 5. Since the central wavelength of the K-band filter (2.20 $`\mu `$m) nearly coincides with the central wavelength of our spectra, no attempt has been made to correct for the different slopes in the continuum. The error introduced by this assumption is well within the accuracy of the flux calibration. The K-band fluxes of four sources have been marked as a lower limit. For IRAS 16279$``$4757 and IRAS 16594$``$4656 the reason is that the absorption profile extends beyond the observed spectral range, hence the estimated continuum is a lower limit to the true level. IRAS 13428$``$6232 and IRAS 17009$``$4154 are larger than the slit at these wavelengths. Flux will have been missed for other sources as well, because the peak-up procedure to position the object in the slit was not very accurate. No attempt has been made to correct for slit loss, because it was impossible to optically verify how well the object was centered in the slit.

When a clear continuum was present, the spectra were normalized by fitting a second order polynomial to the continuum and dividing the spectrum by the fit. The spectra of IRAS 16279$``$4757 and IRAS 16594$``$4656 showed a very wide Br$`\gamma `$ absorption and no or hardly any continuum. Therefore, the spectrum of IRAS 16279$``$4757 was divided by the maximum flux present in the spectrum. For IRAS 16594$``$4656 some continuum seemed present and a linear fit was made between both ends of the spectrum. This normalization should be considered uncertain. The continuum of the four detected PNe is very weak and therefore they were not normalized.

A Voigt profile was used to fit the Br$`\gamma `$ absorption lines in the normalized spectrum. From the fits, the equivalent widths $`W_\lambda `$ were determined using:

$$W_\lambda =_0^{\mathrm{}}[1r(\lambda )]d\lambda ,$$

where $`r(\lambda )`$ is the residual flux normalized to 1 at the continuum. The emission lines were unresolved at our instrumental resolution, and could be fitted well with a Gaussian profile. Table 6 lists the equivalent widths of the absorption and emission profiles, along with the Doppler FWHM $`\mathrm{\Delta }\nu _\mathrm{D}`$, and the Lorentz FWHM $`\mathrm{\Delta }\nu _\mathrm{L}`$ = $`\mathrm{\Gamma }/2\pi `$ (where $`\mathrm{\Gamma }`$ is the effective damping constant) of the Voigt profile.

For IRAS 14488$``$5405 and IRAS 16594$``$4656, both absorption and emission components were present. A combination of a Voigt and a Gaussian profile, each with its own central wavelength, was used in the fit for these objects. For IRAS 16594$``$4656, the central part of the absorption and the emission were fitted well, but the outer wings of the absorption were not. For IRAS 16086$``$5255, a weak absorption line at the central wavelength of Br$`\gamma `$ was observed. It could not fitted well. For IRAS 11159$``$5954, and IRAS 15553$``$5230, no Br$`\gamma `$ was seen in emission or absorption. For IRAS 16279$``$4757 the resulting ‘best fit’ was so bad that no parameters are listed. The fact that the absorption lines in IRAS 16279$``$4757 and IRAS 16594$``$4656 cannot be fitted properly with a Voigt profile, indicates the presence of an additional broadening mechanism. One likely candidate is the linear Stark effect. Further investigations are needed to verify this, but if confirmed, this would indicate a high surface gravity of the central star.

For the spectra containing emission features, the absolute line fluxes were determined. In most spectra only Br$`\gamma `$ was detected, but some also showed the presence of He i emission. In these spectra the the He i $`4f7g`$ $`{}_{}{}^{1}F_{}^{}^1G`$ and $`{}_{}{}^{3}F_{}^{}^3G`$ multiplets would be blended with the Br$`\gamma `$ line, while the He i $`4d7f`$ $`{}_{}{}^{1}D^1F^{}`$ and $`{}_{}{}^{3}D^3F^{}`$ blend could be seen separately. It is unlikely that He ii $`814`$ line was also blended with Br$`\gamma `$, except maybe for high excitation PNe ($`T_{\mathrm{eff}}>\mathrm{80\hspace{0.17em}000}`$ K), such as IRAS 18401$``$1109. In the fitting procedure, it was assumed that each line had a Gaussian profile with the same (unresolved) width and that the wavelength interval between the lines was fixed. The continuum was assumed to be flat except for IRAS 14488$``$5405 and IRAS 16594$``$4656 where the underlying Br$`\gamma `$ absorption was assumed to have a Voigt profile. The helium lines were only fitted when there were indications for their presence upon visual inspection. The separate flux values in the decomposition of the Br$`\gamma `$ complex could not be determined accurately and therefore only the total flux of the blend is listed. The He i lines never contributed more than 8.5 % to the total flux of this blend. The results of these fits are listed in Table 6. No attempt has been made to correct any of the fluxes for slit loss.

The central wavelengths of the hydrogen absorption and emission features can be used to calculate the heliocentric velocities of these objects, using the routine rvcorrect in iraf. The results can be found in Table 6. The literature values for the radial velocities of the standard stars (Hirshfeld et al. Hirsh91 (1991)) were compared with their observed radial velocities, yielding an average accuracy of 11 km s<sup>-1</sup>.

### 3.3 Radio observations

The data were reduced using the package miriad following standard reduction steps as described in the reference guide by Bob Sault and Neil Killeen (http://www.atnf.csiro.au/Software/Analysis/miriad). Images were made using the multi-frequency synthesis technique and robust weighting with a robustness parameter of 0.5. Any confusing sources were clean-ed before determining the upper limits and noise in the map.

We didn’t detect any of the six emission line sources above a detection limit of 0.70 mJy/beam at 3 cm and 0.55 mJy/beam at 6 cm. Our individual upper limits to the flux for the sources are given in Table 7. The weakest confusing source which was clearly detected was 1 mJy and the strongest 3.8 mJy. The map of IRAS 15544$``$5332 showed large scale structure, but still no source was present after deleting the two shortest baselines.

Using the Br$`\gamma `$ flux and assuming Case B recombination we calculated the expected optically thin radio flux. These values are given in the last column of Table 7. The predicted values appear to be a factor ten or more higher than the upper limits on the observed flux. Even the flux per beam of the two extended objects would have been well above the detection limit at 6 cm. For the PNe, however, the predicted values are somewhat lower than what was observed in the radio. The Br$`\gamma `$ flux is probably underestimated due to extinction and slit loss. The radio flux has an uncertainty of 10 % to 20 %.

## 4 Discussion of individual IRAS objects

10256$``$5628: The position of the IRAS counterpart agrees with that of a m(red) = 17.1 mag star in the USNO catalogue. The K magnitude of García-Lario et al. (Garcia97 (1997)) is in agreement with our value.

11159$``$5954: This object has an USNO counterpart with m(red) = 21.0 mag. This is the only star for which the K magnitude derived from the Br$`\gamma `$ spectrum is significantly brighter than the one determined from the K-band image. Hence this star may be variable. Its optical spectrum shows that this is an M-type star (Van de Steene et al., in preparation).

13356$``$6249: Our value for the K magnitude is in perfect agreement with the magnitude mentioned by García-Lario et al. (Garcia97 (1997)). The position and magnitude are also in excellent agreement with the data published by the DENIS project (Epchtein et al. Epch94 (1994)). This object appears extended in our timmi image.

13416$``$6342: According to Hu et al. (Hu93 (1993)), the source is a highly reddened G1I star with R = 17.4 mag. Taking into account the photometry of Hu et al. (Hu93 (1993)), this object seems to have become brighter since the 1987 observations of van der Veen et al. (vdVeen89 (1989)) by $`\mathrm{\Delta }`$J = 0.74 mag, $`\mathrm{\Delta }`$H = 0.46 mag and $`\mathrm{\Delta }`$K = 0.21 mag respectively. However in the L-band the star became fainter by 0.30 mag.

13428$``$6232: This is a very nice bipolar showing Br$`\gamma `$ in emission. The magnitude of van der Veen et al. (vdVeen89 (1989)) is somewhat fainter than ours. Considering its size, possibly part of the flux was missed in their 10″ diaphragm.

13529$``$5934: In the caspir images we see two faint stars less than 4″ apart. The top one is closest to the IRAS position and the reddest source, and therefore the most likely counterpart. The correct identification of the IRAS counterpart needs to be confirmed at 10 $`\mu `$m. The spectra of the objects are completely blended in the irspec image. Their combined spectrum shows no Br$`\gamma `$ in absorption or emission. We won’t discuss this spectrum any further.

14325$``$6428: This object is associated with a bright GSC star of V = 11.8 mag. It was also observed by the DENIS project. Their K<sub>s</sub> magnitude is in good agreement with ours.

14488$``$5405: This IRAS object is also associated with a bright V = 11.5 mag GSC star. This star shows Br$`\gamma `$ in emission.

15066$``$5532: The counterpart at the IRAS position is a bright m(red) = 14 mag USNO star. Its magnitude is K = 9.46 mag. However, our timmi position, and hence the object for which we obtained a spectrum with irspec, is 5$`\stackrel{}{.}`$5 to the west. It is the reddest of the two and its magnitude is K = 8.95 mag. It showed a flat Br$`\gamma `$ spectrum. The object appeared extended in our timmi image. The K magnitude determined by García-Lario et al. (Garcia97 (1997)) is brighter than each one of ours, but he may have had both neighbors in its 15″ aperture. The correct identification of the IRAS counterpart needs to be confirmed at 10 $`\mu `$m. We won’t discuss this spectrum any further.

15144$``$5812: The K magnitude determined in the caspir image is 1 mag brighter than the magnitude determined from the irspec continuum spectrum. We see this bright object and the fainter star towards the south both in the slit. Because of this and the Br$`\gamma `$ in emission, easily visible in the spectrum of the brightest star, there is no doubt that the difference in magnitude is not due to an identification error. Hence this source could be variable.

15553$``$5230: Our K magnitude is 4 mag fainter than the value determined by García-Lario et al. (Garcia97 (1997)). These determinations cannot correspond to the same source. In the K-band image we see four objects close together. The K magnitude of the brightest star to the west is K = 6.21 mag and does agree with the value determined by García-Lario et al. (Garcia97 (1997)). There is a fainter, very red, extended object just eastward of this bright source. It has an elliptical morphology. Most likely this is the correct IRAS counterpart.

16086$``$5255: Our USNO counterpart is significantly fainter than the magnitude K = 7.30 mag determined by Hu et al. (Hu93 (1993)). Their magnitude is in agreement with our magnitude of the bright, red USNO star to the south. Its spectral type would be M3I. Our position is within 3″ of the IRAS position, while the southern star is more than 16″ away. We therefore adopt our identification as the true IRAS counterpart.

16130$``$4620: In the K-band image the object corresponds to two sources 4″ apart. Hu et al. (Hu93 (1993)) associated the southern star (V = 16.7 mag) with an M5Ib star. Our magnitude is in agreement with his measurement. However the top one is by far the reddest of the two, invisible in the optical. This makes it a much stronger candidate for being the IRAS counterpart. The spectra of both objects were blended into one extended source in the irspec observations. We shall not discuss this spectrum any further.

16279$``$4757: Our K magnitude is 0.1 mag fainter and our L magnitude 0.7 mag fainter than the values measured by van der Veen et al. in 1987 and Hu et al. in 1990. However, the J\- and H-band values are in perfect agreement. This difference could be due to measurement errors, or the object may be variable. Hu et al. (Hu93 (1993)) determined R = 18.4 mag and classified it as a G5 star based on its optical spectrum.

16594$``$4656: Our magnitude is in perfect agreement with what was measured by García-Lario et al. (Garcia97 (1997)). HST images show the presence of a bright central star surrounded by a multiple-axis bipolar nebulosity with a complex morphology and a size of 6$`\stackrel{}{.}`$3 $`\times `$ 3$`\stackrel{}{.}`$3 (Hrivnak et al. Hrivnak99 (1999)). In our K-band images we do not see any evidence of the bipolar nebula. The counterpart of this object is a m(red) = 12.4 mag USNO star.

17009$``$4154: Our magnitudes obtained in 1998 are fainter than the 1990 value of García-Lario et al. (Garcia97 (1997)), but brighter than their 1992 values. This source seems variable. It is a post-AGB star with Br$`\gamma `$ in emission surrounded by a faint nebulosity. The object is 7″ in diameter and would have fitted in the aperture of García-Lario, if well centered.

17088$``$4221: This source has been observed by van der Veen et al. (vdVeen89 (1989)) and García-Lario et al. (Garcia97 (1997)). Their K magnitudes are about 1 mag brighter than ours. The position given by van der Veen et al. (vdVeen89 (1989)) corresponds to the bright, red USNO star west of the IRAS position. We measured a magnitude K = 9.08 mag for this star in our image, in agreement with their value. The star for which we obtained the Br$`\gamma `$ spectrum corresponds to the USNO star south of the IRAS position. This object is further away from the IRAS position. The northern star is more likely to be the IRAS counterpart, although the correct identification needs to be confirmed at 10 $`\mu `$m. We won’t discuss this object any further.

## 5 Discussion

### 5.1 Identification of the IRAS counterpart

Due to poor telescope pointing of the ESO 3.6-m at the time of the observations, the counterparts of 3 IRAS sources (13529$``$5934, 15066$``$5532, and 17088$``$4221) could not be determined unambiguously, even after careful analysis of the images. Due to the poor spatial resolution of irspec, the spectrum of IRAS 16130$``$4620 was useless due to blending. We won’t consider these objects any further, so that 16 objects remain for further discussion.

### 5.2 Br$`\gamma `$

Based on the Br$`\gamma `$ spectra we can identify objects of three different types: those with a Br$`\gamma `$ emission line, those with a Br$`\gamma `$ absorption line, and those with no Br$`\gamma `$ line at all. We detected Br$`\gamma `$ in absorption in 7 out of 16 objects. The absorption lines are very narrow in 6 objects, indicating a low surface gravity. This is a strong indication for the post-AGB nature of these objects. Six objects show Br$`\gamma `$ in emission. Two of these also show a photospheric absorption profile. All emission line sources have a strong underlying continuum, unlike normal PNe. In another three objects no clear Br$`\gamma `$ absorption or emission was visible.

### 5.3 Radio continuum

As noted in Sect. 3.3, the predicted optically-thin radio flux values (assuming Case B recombination) appear to be at least a factor of ten higher than the observational flux limits for the post-AGB objects. This indicates that for these objects either the radio flux is optically thick at 3 cm, or the Case B assumption is not valid.

After a star leaves the AGB, its mass loss decreases by several orders of magnitude while simultaneously the velocity of the wind increases. Hence the star will be surrounded by an increasingly more tenuous wind inside a detached AGB shell. At first the star will be too cool to cause any significant ionization, either in the wind or in the AGB shell. However, as the stellar temperature increases, ionization will start and the ionization front will steadily move outwards. This will also be the case for the radius at which the 3 cm radiation becomes optically thick (for brevity we will call this the 3 cm radius, $`r_{3\mathrm{c}\mathrm{m}}`$). Optically thick radio emission from a sphere of radius $`r_{3\mathrm{c}\mathrm{m}}`$ and a temperature T<sub>e</sub> gives a flux density at frequency $`\nu `$ of

$$S_\nu =\frac{r_{3\mathrm{c}\mathrm{m}}^2}{D^2}\frac{2\pi \nu ^2kT_\mathrm{e}}{c^2}.$$

where D is the distance to the star.

Using this approximation we calculated that $`r_{3\mathrm{c}\mathrm{m}}<`$ 3$`\times 10^{14}`$ cm, assuming $`S_{3\mathrm{c}\mathrm{m}}<`$ 0.7 mJy, $`T_\mathrm{e}`$ = $`10^4`$ K, and $`D`$ = 1 kpc. This upper limit is similar to the size Knapp et al. (Knapp95 (1995)) determined for the post-AGB stars CRL 915 ($`S_{3\mathrm{c}\mathrm{m}}=`$0.30 mJy) and IRAS 17423$``$1755 ($`S_{3\mathrm{c}\mathrm{m}}=`$0.44 mJy) using the same assumptions. As long as the ionization front has not reached the AGB shell, the 3 cm radius will not change much, because the outer regions in the $`1/r^2`$ density profile contribute very little to the optical depth, and consequently the radio flux will remain very low. However, once the post-AGB star reaches a temperature hot enough to ionize the AGB shell, the 3 cm radius will quickly increase by roughly two orders of magnitude, causing a dramatic increase in radio flux. This marks the onset of the PN phase. As an example, the sizes of the young PNe CRL 618 and IRAS 21282+5050 measured by Knapp et al. (Knapp95 (1995)) are a factor 10 larger than the sizes of post-AGB star nebulae. Their radio flux values, 67 mJy and 4.3 mJy respectively, are well above our detection limit.

In our post-AGB star candidates with Br$`\gamma `$ in emission, the ionized region where the emission originates must be very small and dense. Probably, the AGB shell of our objects is not yet ionized, but the post-AGB wind could be. The Br$`\gamma `$ spectra are unresolved at a resolution of $``$150 km s<sup>-1</sup>. Hence the wind velocity couldn’t be much larger than this value. Evidence for a wind emanating from some of these central stars was presented in Van de Steene et al. (VdSteene00 (2000)).

### 5.4 Spectral energy distribution

One of the well-established characteristics of post-AGB stars is that their Spectral Energy Distributions (SEDs) have a ‘double-peaked’ shape. The two peaks in the spectrum correspond to the stellar and dust emission components. Post-AGB stars have been classified into four classes based on the shape of the SED by van der Veen et al. (vdVeen89 (1989)).

* class I: has a flat spectrum between 4 $`\mu `$m and 25 $`\mu `$m and a steep fall-off to shorter wavelengths.
* class II: maximum around 25 $`\mu `$m and a gradual fall-off to shorter wavelengths.
* class III: maximum around 25 $`\mu `$m and a steep fall-off to a plateau roughly between 1 $`\mu `$m and 4 $`\mu `$m with a steep fall-off at shorter wavelengths.
* class IV: two distinct maxima, one around 25 $`\mu `$m and a second between 2 $`\mu `$m and 3 $`\mu `$m (IVa), between 1 $`\mu `$m and 2 $`\mu `$m (IVa), or below 1 $`\mu `$m (IVb).

Note that class IVa was not contained in the original definition, but was added to classify objects that did not fit in any of the original categories.

The SEDs of the objects are shown in Fig. 10. The objects have the typical post-AGB SEDs as cited above. The SED class of each object is listed in Table 5. Six of the 16 positively identified objects are of class II and the other 10 of class IV. For objects in class II the circumstellar dust is so optically thick that almost all star light is absorbed by the dust and is re-radiated at mid- to far-infrared wavelengths. The large infrared excess is commonly attributed to the presence of a very compact circumstellar dust shell and/or ongoing mass loss which obscures the central star from view. Objects in class IV have less obscured central stars: the thermal emission from their circumstellar shells appear as a peak in the far-infrared and the central stars show up as a peak in the near-infrared (IVa) or optical (IVb). We always see some stellar signature in the near-infrared and therefore have no objects of class III. For instance, Van der Veen et al. (vdVeen89 (1989)) classified IRAS 16594-4656 as class III, while we classified it as IVa. We extended the definition of class IVa to include objects which show a clear stellar signature in the near-infrared, but peak in the K-band, just beyond 2 $`\mu `$m. (e.g. IRAS 13428-6232, IRAS 16279-4757). In the table these objects are marked as class IVa. These objects have no optical counterpart in the USNO catalog. When objects in class IV peak at shorter wavelengths they usually have an USNO counterpart. The central stars of objects in class IVb are less reddened and brighter than objects in class IVa and have bright optical counterparts (e.g. IRAS 14325-6428, IRAS 14488$``$5405).

We especially draw attention to two objects which have unique SEDs. IRAS 15544$``$5332 is the only object in the sample for which the L-band value is higher than the IRAS 12 $`\mu `$m value. It has the coolest dust shell in the sample. It also has a very steep class II spectrum, indicative of a very high extinction. IRAS 11159-5954 is the only object in the sample for which the peak of its SED in the near-infrared is higher than in the far infrared, showing that the grain emission is very weak.

### 5.5 Color-color diagrams

#### 5.5.1 IRAS color-color diagram

In Fig. 1 we show the IRAS color-color diagram. The IRAS fluxes were converted to magnitudes according to the IRAS Explanatory Supplement (Beichman et al. Beichman84 (1984)). The boxes defined by van der Veen & Habing (vdVeen88 (1988)) are drawn in. According to this classification scheme, PNe are found in region V of the color-color diagram and AGB stars in region IV. In region VIII there may be some confusion from galaxies and young stellar objects, and in region IV an odd H ii region may be present. The only object in region VIII is IRAS 15544$``$5332, which is not redshifted and shows Br$`\gamma `$ in emission. The object is unresolved, and the Br$`\gamma `$ emission is very weak. Hence it is unlikely to be an ultra-compact H ii region, but we cannot completely rule out that it is an embedded young stellar object that is not hot enough to ionize its surroundings. IRAS 13416$``$6243 in region IV has Br$`\gamma `$ in absorption and hence can be considered to be a post-AGB star. ¿From the results in Fig. 1, it seems that in the IRAS color-color diagram no distinction can be made between post-AGB stars with Br$`\gamma `$ in emission, absorption, or a flat Br$`\gamma `$ spectrum. Post-AGB stars were expected to be located in a region in the IRAS color-color diagram between AGB stars and PNe (e.g. Volk & Kwok Volk89 (1989); van der Veen et al. vdVeen89 (1989); Hu et al. Hu93 (1993)). However van Hoof et al. (vHoof97 (1997)) found in their parameter study of the spectral evolution of post-AGB stars that they can follow a variety of paths in the IRAS color-color diagram. Consequently PNe and post-AGB stars can occupy the same region in the IRAS color-color diagram and the position in the IRAS color-color diagram alone cannot a priori give a unique determination of the evolutionary status of a post-AGB star. Our observations confirm this result.

The fact that our objects are mostly selected from the region where PNe were found and not in the region between AGB and PNe, may explain our high detection rate of Br$`\gamma `$ emission (for comparison, the search by Käufl et al. (Kaufl93 (1993)) for Br$`\alpha `$ emission in a sample of 21 post-AGB stars resulted in only one detection).

#### 5.5.2 Near-infrared color-color diagrams

In Fig. 2 and Fig. 3 we show the J$``$H versus H$``$K and the H$``$K versus K$``$L color-color diagrams, respectively. In the diagrams we notice that objects of class II are redder than the objects of class IV. The former are found in the top right part of the diagrams while the latter are located more towards the bottom left. The separation is most obvious in Fig. 3. For our sample, objects having K$``$L $`>`$ 1.5 mag are all of class II and objects which have a maximum in the near-infrared beyond 2 $`\mu `$m have H$``$K $`>`$ 1.0 mag. Objects for which the stellar signature is more pronounced and which peak more towards shorter wavelengths in the near-infrared, are found more towards the bottom left in the diagrams.

To understand this effect, we would like to point out that for all but the coolest stars the intrinsic shape of the stellar continuum in the near-infrared can be approximated by a Rayleigh-Jeans tail. Since the shape of this tail does not depend on stellar temperature, the intrinsic infrared colors of these stars will not depend on stellar temperature either. Furthermore, since the colors of an A0V star are by definition zero, the infrared colors for most other stars are also close to zero, except when severe inter- or circum-stellar extinction is present. In this case the infrared colors will be positive and can be used as a crude measure for the extinction.

The mass loss rate during the superwind phase is very high. At the end of this phase the circumstellar dust will almost completely obscure the central star. During the post-AGB evolution, as the AGB shell expands, the shell will become more optically thin to stellar radiation. The stellar signature becomes more pronounced and will peak towards ever shorter wavelengths in the near-infrared. Consequently, it is expected that post-AGB stars will move from the upper right of the diagram to the lower left as the AGB shell expands and the circumstellar extinction becomes less. Generally speaking, we might expect that the objects in the top right of the diagram have left the AGB more recently than the objects in the lower left. However, we need to be cautious about such an interpretation because of the combined effects of many unknowns. The time it takes for the envelope to become optically thin will depend upon the mass loss rate at the tip of the AGB, the wind velocity, and the distribution of the mass in the circumstellar shell. If the AGB star had a non-spherical mass loss concentrated towards the equator, the central star could be visible along polar directions, while being completely obscured in equatorial directions. Hence, the observed amount of circumstellar extinction will depend on the viewing angle (Soker Soker99 (1999)). The interstellar extinction also affects the position of the objects in the color-color diagrams, as indicated by the arrow. The J$``$H versus H$``$K diagram is clearly the most affected by interstellar extinction. Eleven of the 16 objects are within one degree of the galactic plane, including the 6 objects in class II. The two objects in class IVb also have the highest latitude ($`|b|>`$ $`4^{}`$). We calculated the extinction for our objects at a distance of 1 kpc and 4 kpc according to Hakkila et al. (Hakkila97 (1997)): $`A_V`$ would be between 0.5 mag and 2.0 mag if the objects were at a distance of 1 kpc and between 2.1 mag and 6.2 mag with a median of 4.6 mag if they were at 4 kpc.

Objects with Br$`\gamma `$ in emission, absorption, or a flat spectrum are well mixed in the diagrams. As discussed in the previous section it is unlikely that the stars with Br$`\gamma `$ in emission have reached a temperature high enough to start to ionize the AGB shell and the emission probably originates in the stellar wind. Because we observe Br$`\gamma `$ in emission from objects in class II (e.g. IRAS 15144$``$5812), it seems that a fast stellar wind can be present at an early stage when the circumstellar shell is still very optically thick. We are currently trying to determine spectral types of our objects in order to determine their true post-AGB evolutionary status (Van de Steene et al., in preparation).

#### 5.5.3 Combined near- and far-infrared color-color diagrams

In Fig. 4 and Fig. 5 we show the K$``$L versus $``$ and K$``$L versus $``$ diagrams, respectively. The K$``$L color roughly describes the evolution of the circumstellar extinction and is therefore a measure for the expansion of the AGB shell, while the $``$ and $``$ colors reflect the evolution of the spectrum of the circumstellar dust shell.

In Fig. 4 we see a weak trend towards cooler $``$ colors with decreasing K$``$L color at first. When K$``$L $`<`$ 1.5 mag this trend seems to reverse, probably due to an increase in 12 $`\mu `$m flux. However, this trend will need to be confirmed with a larger sample.

In Fig. 5 we see a broad and weak trend towards hotter $``$ colors with decreasing K$``$L color. Obviously the 25-$`\mu `$m flux increases faster than the 60-$`\mu `$m flux.

Theoretical calculations by Blöcker (Bloecker95 (1995)) predict that directly after the star leaves the AGB, the stellar evolution is slow. This causes the AGB dust shell to cool when it expands. However, around $`T_{\mathrm{eff}}8000`$ K the evolution of the central star speeds up considerably and the grains in the AGB shell start to heat up again. This effect is most pronounced for silicate grains because, as the peak of the stellar spectrum moves into the UV, the efficiency with which these grains absorb light increases significantly. This causes the counter-clockwise loop which was first predicted by van Hoof et al. (vHoof97 (1997)) for oxygen-rich post-AGB stars in the IRAS color-color diagram. This effect could explain the reverse trend in Fig. 4, and also the trend in Fig. 5.

In Fig. 6 we plot the K$``$L versus L$``$ diagram. The K$``$L color describes the evolution of the circumstellar extinction. The L$``$ color relates the stellar component with the peak of the dust emission.

Because he K$``$L color is sensitive to the extinction and the L$``$ color is insensitive to extinction, all class II objects are found in the upper half of the diagram and class IV objects in the lower half.

For our limited sample, all but one of the objects are found to the right of the dotted line. This is partly due to sample selection. As the arrow indicates, objects less obscured (e.g., due to orientation along the polar axis, or because of an optically thin shell, or less interstellar extinction), and especially the ones with cool central stars, could be situated below the dotted line.

IRAS 11159$``$5954 has a lower extinction than its L$``$ color would indicate. This is an M-type star (Van de Steene et al., 2000, in preparation), which may still have ongoing mass loss. Its dust shell does not appear to be very thick and the star is very bright in the near infrared.

The objects to the right of the dashed line are extended in the near-infrared or the optical. IRAS 17009$``$4154 and IRAS 15553$``$5230 showed elliptical morphology. The former is very faint in the optical, the latter invisible. IRAS 13428$``$6232 shows a bipolar morphology in the near-infrared and and is also very faint in the optical. IRAS 16594$``$4656 is bright and was not observed to be extended in the near-infrared, though it showed a bipolar morphology in its HST image (Hrivnak et al. Hrivnak99 (1999)). Probably they have a thicker cool circumstellar dust shell than objects located below the dashed line (more 25 $`\mu `$m-band flux) and/or their central star temperatures are higher (less L-band flux).

As the dust shell expands, the  magnitude will decrease a bit (for oxygen-rich post-AGB stars) or remain roughly constant (for carbon-rich post-AGB stars) (van Hoof et al. vHoof97 (1997)). The L magnitude will increase, as the stellar temperature increases. The K magnitude will decrease as the star starts shining through the dust shell. Thus the L$``$ values are expected to increase with decreasing K$``$L values as the shell expands and the star becomes hotter.

This is what we observe in Fig. 6, both for the extended and non-extended objects. A small L$``$ color indicates that the system is young, hence the star is heavily obscured by the dust shell which is still close to the star. This is noticeable by the large K$``$L values. For objects with a larger L$``$ color the AGB shell is already more detached and cooler. The extinction is less and this translates into a smaller value for the K$``$L color.

The K$``$L versus L$``$ diagram (not shown) is very similar to the K$``$L versus L$``$ diagram, but the spread is larger for K$``$L $`<`$ 1.0 mag, possibly because of an increase in 12-$`\mu `$m flux, as discussed in the previous section. Because of its relatively large 12-$`\mu `$m flux, IRAS 13416$``$6342 is located close to IRAS 17009$``$4154. It needs to be checked at higher resolution whether this source is extended.

In Fig. 7 we show the color-color diagram proposed by Ueta et al. (Ueta00 (2000)). Since the J$``$K color measures the circumstellar extinction, in a similar way as K$``$L color does, and because the evolution described above for the L$``$ color is equally valid for the K$``$ color, one would expect the same decreasing trend as in the previous diagram. However, extinction effects have a much stronger effect on the K$``$ and J$``$K colors than on the L$``$ and K$``$L colors, as can be seen from the arrows in both diagrams. Nearly all objects of class IVa have J$``$K $`<`$ 2.5 mag. When the star becomes prominent in the near-infrared, the stellar maximum shifts through the K-band towards the J-band, causing a reversal in the evolutionary trend.
The extended objects are also found towards the right of the diagram. In this plot it is more obvious that IRAS 16594$``$4656 is the brightest of the 4 extended objects.

Ueta et al. (Ueta00 (2000)) had one M star, IRAS 04386$`+`$5722, offset towards blue K$``$ color, similar to IRAS 11159$``$5954, but a bit bluer in J$``$K. Its position is indicated with the smallest grey box in Fig. 7.

The long-dashed line at J$``$K = 1.45 mag in Fig. 7 separates the regions what Ueta et al. call Star-Obvious Low level Elongated (SOLE) and DUst Prominent Longitudinally EXtended (DUPLEX) nebulae. One quarter of our objects would be classified as SOLE and three quarters as DUPLEX in this scheme. The grey regions show where the objects in their sample are located and are labeled with their acronyms. Our samples are obviously complementary: there is virtually no overlap ! The difference may be caused in part by the selection criteria. They have selected known post-AGB candidates from the literature and imaged their nebulosities with WFPC2 in the optical. Consequently they chose optically bright post-AGB stars. We selected objects from the PNe region in the IRAS color-color diagram of which very few had an optical counterpart identified. Moreover, all our objects are located within 5 degrees of the galactic plane, while only 11 out of their 27 objects are (3 SOLE, 6 DUPLEX, and 2 stellar). The arrow in the diagram indicates a correction for $`A_V`$ = 5 mag of extinction. Larger interstellar extinction alone cannot explain why our samples appear different. If the objects in both samples are at similar distances, it is plausible that we have more massive central stars in our sample, and that they are evolving faster across the HR diagram. Further investigation is needed to understand the differences between both samples.

In summary, no distinction can be made between the objects showing Br$`\gamma `$ in emission, absorption, or a flat spectrum, in any of the color-color diagrams. The trends we see in the near and far infrared are mainly due to the expansion, morphology, and dust properties in the circumstellar shell and the obscuration of the central star it causes. The trends show the expected evolution of the circumstellar shell. Whether the positions of the objects in the color-color diagrams can be directly related to the temperature and core mass of the central star needs further investigation.

## 6 Conclusions

In this article we reported further investigations of the IRAS selected sample of PN candidates that was presented in Van de Steene & Pottasch (VdSteene93 (1993)). About 20 % of the candidates in that sample have been detected in the radio and/or H$`\alpha `$ and were later confirmed as PNe. Here we investigated the nature of the non-radio-detected sources.

* Of sixteen positively identified objects, seven show Br$`\gamma `$ in absorption. The absorption lines are very narrow in six objects, indicating a low surface gravity. This is a strong indication for the post-AGB nature of these objects. Another six objects show Br$`\gamma `$ in emission. Two of these also show photospheric absorption lines. All emission line sources have a strong underlying continuum, unlike normal PNe. In another three objects, no clear Br$`\gamma `$ absorption or emission was visible.
* The objects showing Br$`\gamma `$ in emission were re-observed in the radio continuum with the Australia Telescope Compact Array. None of them were detected above a detection limit of 0.55 mJy/beam at 6 cm and 0.7 mJy/beam at 3 cm, while they should have been easily seen if the radio emission was optically thin and Case B recombination was applicable. It is argued that the Br$`\gamma `$ emission may be due to ionization in the post-AGB wind, present before the star is hot enough to ionize the AGB shell.
* The fact that our objects were mostly selected from the region in the IRAS color-color diagram where typically PNe are found, may explain our higher detection rate of emission line objects compared to previous studies, which selected their candidates from a region between AGB and PNe. These post-AGB stars also cover a larger range in color and are generally much redder than the ones known so far.
* In the near-infrared color-color diagrams our objects cover a very large range of extinction. Near-infrared versus far-infrared color-color diagrams show trends which reflect the expected evolution of the expanding circumstellar shell. No distinction can be made between the objects showing Br$`\gamma `$ in emission, absorption, or a flat spectrum in the near- and far-infrared color-color diagrams. Whether the positions of the objects in the color-color diagrams can be directly related to the temperature and core mass of the central star needs further investigation.
* We identified the K$``$L versus L$``$ diagram as a potentially useful tool to distinguish: 1) extended from unresolved post-AGB stars, and 2) obscured objects of class II having thick circumstellar shells from the brighter class IV objects which show a stellar signature in their near-infrared SEDs. However, this result should be confirmed with a larger sample.

###### Acknowledgements.

PvH wishes to thank ESO for their hospitality and financial support during his stay in Santiago where part of this paper was written. PvH acknowledges support by the Netherlands Foundation for Research in Astronomy (ASTRON) through grant no. 782–372–033 during his stay in Groningen and the NSF through grant no. AST 96–17083 during his stay in Lexington. SkyCat was developed by ESO’s Data Management and Very Large Telescope (VLT) Project divisions with contributions from the Canadian Astronomical Data Center (CADC).