A&A 376, 165-174 (2001)
DOI: 10.1051/0004-6361:20010938

The residual orbital eccentricity in close binaries

An application to millisecond binary pulsars

A. F. Lanza1 - M. Rodonò1,2

1 - Osservatorio Astrofisico di Catania, Via Santa Sofia, 78, Città Universitaria, 95125 Catania, Italy
2 - Dipartimento di Fisica e Astronomia dell'Università degli Studi, Via Santa Sofia, 78, Città Universitaria,
95125 Catania, Italy

Received 7 February 2001 / Accepted 27 June 2001

Density fluctuations in convective stars produce small variations of their outer gravitational field that give rise to a very small, though finite, orbital eccentricity in close binaries with convective components. We revisit the theory for such an effect, originally presented by Phinney (1992), and propose a new approach to compute the fluctuations of the outer gravitational field based on Chandrasekhar's virial theorem. It is applied to compute the eccentricities of binaries with a millisecond pulsar component in the framework of a simplified model of their evolution. The comparison with the observations allows us to test the models proposed for the formation of such systems, giving support to the evolutionary scenarios that connect them with low-mass X-ray binaries, and to provide a new test for turbulent convection models of giant stars.

Key words: stars: close binaries - binary pulsars - celestial mechanics - stellar convection - extrasolar planets

1 Introduction

Tidal interactions are particularly effective in close binaries with convective components and lead to orbit circularization and spin-orbit synchronization over time scales much shorter than the evolutionary time scales of the components. However, perfect circularization and synchronization is never attained, as noticed by Phinney (1992). He pointed out that convection produces random fluctuations of the density inside a star, associated with rising and falling turbulent plasma columns, that induce, in turn, fluctuations of the star's outer gravitational field. The random radial force acting between the components excites small oscillations of their separation with the period of the unperturbed orbit. The resulting orbital eccentricity grows proportionally to t1/2, where t is time, in a random walk process, analogue to the resonant excitation of a linear oscillator by a random force. However, its indefinite increase is prevented by tidal dissipation, which grows linearly with the eccentricity itself, until a statistical equilibrium is reached between the stochastic excitation by the convective fluctuations and the tidal dumping.

Phinney noticed that the measurement of the small residual eccentricity produced by this process gives information on the response of the orbit to the tidal dissipation, thus introducing an analogue of the fluctuation-dissipation theorem in statistical physics (cf., e.g., Reichl 1980). The most interesting application is the prediction of the final eccentricities of binaries containing millisecond pulsars, which can actually be measured thanks to the very accurate timing that is now achievable for such systems (cf. Phinney & Kulkarni 1994). The comparison between theory and observations thus provides a powerful test of the evolutionary scenarios proposed for the formation of these binaries.

In the present paper we revisit Phinney's suggestion that convective motions can induce fluctuations of the gravitational field of a convective star and show how, by applying Chandrasekhar's virial theorem, the power spectrum of such fluctuations can be computed. The effects of these fluctuations on the orbital dynamics are evaluated using the methods of general perturbation theory. Only the fluctuations of the quadrupole term in the multipole expansion of the gravitational potential are relevant because higher order multipoles decay very rapidly with the distance from the baricentre of the star. An illustrative application of our approach supports Phinney's conjecture and shows how the distribution of the final eccentricities depends on the assumed properties of convection and tidal dissipation efficiency. Conversely, this provides a novel method to test stellar convection and tidal dissipation.

2 Fluctuations of the gravitational quadrupole moment

Phinney (1992) showed that the gravitational quadrupole moment of a star with a convective envelope is not constant, being subject to fluctuations related to the density fluctuations in the convective envelope. He estimated the amplitude and the power spectrum of the fluctuations of the quadrupole moment in a simplified way. In this Section, we present a more rigorous treatment based on an application of the tensor virial theorem in magnetohydrodynamics after Chandrasekhar (1961). The main advantage of the present approach is that the fluctuations of the quadrupole moment tensor are directly related to the kinetic energy fluctuations produced by convective motions.

In order to study the fluctuations due to convection, we shall consider an inertial reference frame with the origin at the baricentre of the star, the z axis along its rotation axis and the x and y axes lying in the equatorial plane.

The external gravitational field produced by the convective component can be expanded up to the quadrupole term as:

 \begin{displaymath}\Phi({\vec x}) = \frac{G M}{r} + \frac{3}{2} G Q_{ik} \frac{x_{i}x_{k}}{r^{5}}
\end{displaymath} (1)

where ${\vec x}$ is the position vector, r the radial distance from the baricentre of the star, M its mass, G the gravitational constant, Qik is the quadrupole moment tensor, which is defined as the traceless part of the inertia tensor Iik:

 \begin{displaymath}Q_{ik} = I_{ik} - \frac{1}{3} \delta_{ik} {\rm Tr} I,
\end{displaymath} (2)

with the inertia tensor defined as

 \begin{displaymath}I_{ik} = \int_{V} \rho x_{i} x_{k} {\rm d} {\vec x},
\end{displaymath} (3)

where $\rho$ is the density and V the volume of the star (cf. Chandrasekhar 1961).

The centrifugal and tidal distortions induce a static quadrupole moment, which can be neglected for the present purposes although it may be slowly modulated over time scales of several decades in stars with intense magnetic activity (cf., e.g., Applegate 1992; Lanza & Rodonò 1999).

The tensorial virial theorem in our inertial reference frame is expressed by the relationship (Chandrasekhar 1961):

 \begin{displaymath}\frac{1}{2} \ddot{I}_{ik} = 2 {\cal T}_{ik} - 2 {\cal M}_{ik} +{\cal G}_{ik} + \delta_{ik} \{ (\gamma-1){\cal U} + {\cal M} \}
\end{displaymath} (4)

where the inertia tensor Iik is defined by Eq. (3), and

 \begin{displaymath}{\cal T}_{ik} = \frac{1}{2} \int_{V} \rho v_{i} v_{k} {\rm d}{\vec x},
\end{displaymath} (5)

\begin{displaymath}{\cal M}_{ik} = \frac{1}{8\pi} \int_{V} B_{i} B_{k} {\rm d}{\vec x},
\end{displaymath} (6)

$\displaystyle {\cal G}_{ik} =
-\frac{1}{2} G \int_{V}\int_{V}
\rho({\vec x}) \r...
...\vec x} - {\vec x^{\prime} }\vert^{3}} {\rm d}{\vec x} {\rm d}{\vec x^{\prime}}$     (7)

are the kinetic energy tensor, the magnetic energy tensor and the gravitational energy tensor, respectively. In Eqs. (5)-(7), ${\vec v}$ is the velocity of the plasma, ${\vec B}$ the magnetic field, V the volume of the star, $\gamma$ the ratio of the specific heats, G the universal constant of gravitation and ${\cal U}$ the internal energy of the star (the time derivative is indicated by a dot).

From the definition of the gravitational quadrupole moment (Eq. (2)), Eq. (4) implies:

$\displaystyle \ddot{Q}_{ik}$=$\displaystyle 4\left[ \left( {\cal T}_{ik} - \frac{1}{3}\delta_{ik}{\cal T}
\right) - \left( {\cal M}_{ik} + \frac{1}{6}\delta_{ik}{\cal M}\right) \right]$
$\displaystyle + 2 \left( {\cal G}_{ik} -\frac{1}{3} \delta_{ik}{\cal G} \right),$ (8)

where the terms without indexes indicate the traces of the corresponding tensors. It is useful to group together the terms containing the gravitational energy tensor by defining:

 \begin{displaymath}{\Xi}_{ik} \equiv {\cal G}_{ik} - \frac{1}{3} \delta_{ik} {\cal G}.
\end{displaymath} (9)

Hereinafter we shall be concerned with random fluctuations of the quantities that appear in the above equations and shall denote such fluctuations with a prime. In Appendix A it is shown that the fluctuation of $\Xi_{ik}$, under quite general hypotheses, can be written as: $\Xi_{ik}^{\prime}
= -\omega_{\rm f}^{2} Q_{ik}^{\prime}$, where $\omega_{\rm f}$ is a characteristic frequency that depends on the amplitude of the density fluctuations and the unperturbed density distribution inside the star. The equation for the fluctuations of the quadrupole moment tensor can then be derived from Eq. (8):

 \begin{displaymath}\ddot{Q}^{\prime}_{ik} + \omega_{\rm f}^{2} Q ^{\prime}_{ik} ...
...} + \frac{\delta_{ik}}{6}{\cal M}^{\prime}\right) \right]\cdot
\end{displaymath} (10)

Equation (10) was obtained from the virial theorem of ideal magnetohydrodynamics which does not include dissipative processes. In a real star, dissipative processes damp the fluctuations of the quadrupole moment and we may take them into account by adding a damping term to the l.h.s. of Eq. (10) (cf. also Chandrasekhar 1969). The characteristic damping time $\tau_{\rm d}$ will be of the order of the turbulent dissipation time $\tau_{\rm d} \sim R^{2}/\nu_{\rm t}$, where R is the star's radius and $\nu_{\rm t}$ the turbulent viscosity. Note that in a convection zone $\tau_{\rm d}$ is pratically the same for the dissipation of a shear flow and for the decay of a magnetic field because the turbulent coefficients for diffusion processes are all of the order $\nu_{\rm t} \sim \ell_{0} u_{0}$, where $\ell_{0}$ is the mixing length and u0 the turbulent velocity of the largest convective eddies. Therefore, the equation governing the fluctuations of the quadrupole moment can be written as:

\begin{displaymath}\ddot{Q}^{\prime}_{ik} + \lambda \dot{Q}^{\prime}_{ik} + \omega_{\rm f}^{2} Q^{\prime} _{ik} = \end{displaymath}

$\displaystyle \ \ \ \ \ \ 4\left[ \left( {\cal T}^{\prime}_{ik} - \frac{\delta_...
... {\cal M}^{\prime}_{ik} + \frac{\delta_{ik}}{6}{\cal M}^{\prime}\right) \right]$     (11)

where $\lambda = 1/\tau_{\rm d}$.

In order to study the fluctuations of the energy terms, we express the instantaneous density and velocity at a given position ${\vec x}$ inside the star as:

 \begin{displaymath}\rho ({\vec x}, t) = \rho_{0} ({\vec x}) + \rho^{\prime} ({\vec x}, {\it t})
\end{displaymath} (12)

 \begin{displaymath}{\vec v({\vec x}, {\it t})} = {\vec v}_{0}({\vec x}) + {\vec v}^{\prime}({\vec x}, t)
\end{displaymath} (13)

where $\rho_{0}$ and ${\vec v}_{0}$ are the mean density and velocity, respectively, and $\rho^{\prime}$, ${\vec v}^{\prime}$ are their fluctuations, whose time and ensemble averages are assumed to be zero. We shall also assume that the ensemble averages of the products of the velocity components, i.e. their correlations, are given by:

 \begin{displaymath}\langle v^{\prime}_{i} v^{\prime}_{k} \rangle
= \delta_{ik} \langle v^{\prime \, 2} \rangle .
\end{displaymath} (14)

This implies that the correlation between different velocity components vanishes. As a matter of fact, in a rotating star there may be some degree of correlation between different components of the velocity, but the effect can be neglected for our purpose. Substituting Eqs. (12) and (13) into the definition of the kinetic energy tensor (Eq. (5)) and performing the integration over the volume of the star, the first order fluctuating components are:

 \begin{displaymath}{\cal T}_{ik}^{\prime} = \frac{1}{6}\delta_{ik} \int_{V} \rho_{0} v^{\prime \,2} {\rm d}{\vec x}
\end{displaymath} (15)

where the hypotheses on the correlations between the fluctuating terms have been applied. The fluctuations of the magnetic energy may be considered much smaller than those of the kinetic energy and will be neglected in Eq. (11).

In view of Eq. (15), the only non-trivial components in Eq. (11) are those with i=k. The fluctuations of the gravitation quadrupole moment tensor can be characterized by means of its power spectrum PQik (see, e.g., Panter 1965 for a definition of the power spectrum and its basic properties). Substituting Eq. (15) into Eq. (11) and considering that the statistical properties of the fluctuations are spherically symmetric, we find:

 \begin{displaymath}P_{Q^{\prime}_{xx}} (\omega) = \frac{32}{9} \frac{P_{{\cal T}...
...\omega_{\rm f}^{2} -
\omega^{2})^{2} + \lambda^{2}\omega^{2}},
\end{displaymath} (16)

where $\omega$ is the frequency and $P_{{\cal T}}(\omega)$ is the power spectrum of the kinetic energy fluctuations. Similar equations hold for the power spectra of the other non-vanishing components $Q_{yy}^{\prime}$ and $Q_{zz}^{\prime}$.

The determination of $P_{{\cal T}}(\omega)$ calls for a model of the convective turbulence in stellar interiors. The Reynolds and the Rayleigh numbers in stellar convection are many orders of magnitude larger than in laboratory experiments or in detailed numerical simulations, thus significant uncertainties remain on the statistical properties of stellar convective flows (cf., e.g., Schatzman & Praderie 1993; Rüdiger 1989). Nevertheless, we shall assume a simplified model for stellar convection in order to estimate the power spectrum of the kinetic energy fluctuations. For simplicity sake, we subdivide the convective envelope into domains whose fluctuations are statistically uncorrelated; as a consequence the power spectrum of the entire convection zone is simply the sum of the power spectra of the single domains (cf. Panter 1965). Moreover, we shall assume that the density is constant within a given domain, which is a fairly good assumption since the typical length scale of the largest convective eddies is assumed to be comparable with the local pressure (and density) scale height, according to the standard mixing-length treatment.

In convection models based on the mixing length approach, the energy is injected in the largest eddies and then it goes down along a turbulent cascade, until the smallest eddies reach the dimension at which viscosity effects become important. The energy cascade produces a spectrum which is of the Kolmogorov's type (e.g., Goldreich & Keeley 1977; Goldreich & Kumar 1988). If we denote $\ell_{0}$ and u0 as the length scale and the velocity of the largest eddies, the velocity $v(\omega)$ of an eddy of turnover frequency $\omega$ is $v(\omega) \sim u_{0} (\omega/\omega_{0})^{-1/2}$, where $\omega_{0} = u_{0}/\ell_{0}$. The ratio of the mixing length $\ell (\omega)$ of that eddy to the length of the largest eddies $\ell_{0}$ scales as: $\ell (\omega)/ \ell_{0} \sim (\omega/\omega_{0})^{-3/2}$ and its kinetic energy $ e(\omega)$ - not to be confused with the power spectrum of the kinetic energy - as $ e(\omega) \sim e_{0} (\omega/\omega_{0})^{-11/2}$, where e0 is the kinetic energy of the largest eddies. The fluctuations of the kinetic energy of the largest eddies of dimension $\ell_{0}$due to the eddies of dimension $\ell (\omega)$ scale as  ${\sim} (\omega/\omega_{0})^{-13/2}$ because the largest eddies contain $(\ell_{0}/\ell)^{3}$ uncorrelated smaller eddies of dimension $\ell (\omega)$ and energy $ e(\omega)$. Therefore, the power spectrum of the kinetic energy fluctuations in the inertial range, where the Kolmogorov spectrum applies, scales as $P_{\cal T}(\omega)
\sim (\omega/\omega_{0})^{-15/2}$.

The Kolmogorov spectrum can not be extended to very low frequencies ( $\omega \ll \omega_{0})$ because that would imply an infinite total energy per unit volume. Therefore, for $\omega < \omega_{0}$we shall assume that the fluctuations of the total kinetic energy can be regarded as hydrodynamics fluctuations (cf., e.g., Lifshitz & Pitaevskii 1980), with a characteristic relaxation time $\tau_{\rm r}$ equal to the turbulent dissipation time of the layer: $\tau_{\rm r} =
\ell_{0}^{2}/\nu_{\rm t} \sim \ell_{0}^{2}/(u_{0}\ell_{0}) \sim \omega_{0}^{-1}$.

In the light of these hypotheses, the power spectrum of the kinetic energy fluctuations can be written as:

 \begin{displaymath}P_{\cal T} (\omega) = \left\{
\frac{ 4 \o...
... & \mbox{ for $\omega \geq \omega_{0}$ }\\
\end{displaymath} (17)


\begin{displaymath}P_{\cal T} (\omega) = P_{\cal T} (-\omega) \;\;\;\;\;\;\;\;\;...
...$\omega < 0$ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

because ${\cal T}^{\prime}$ is real. It has been normalized in such a way that: $\int_{-\infty}^{\infty} P_{\cal T} (\omega) {\rm d}\omega \simeq 2
... P_{\cal T} (\omega) {\rm d}\omega = 2 \pi \langle {\cal T}^{\prime\,2}\rangle $. The mean square fluctuation of the kinetic energy of a layer of radius r and mixing length $\ell$ is $\langle {\cal T}^{\prime\,2}\rangle \sim 4\pi (r/\ell)^{2} {e}^{2}_{0}$.

The power spectrum of the quadrupole moment fluctuations produced by a layer of mixing length $\ell$ can be easily found by substituting Eq. (17) into Eq. (16). The power spectrum of the quadrupole moment fluctuations of the whole star is obtained by summing all the contributions coming from the uncorrelated convective layers, into which the convection zone was divided.

The mean square fluctuation of the quadrupole moment produced by a convective layer can be computed by integrating Eq. (16):

 \begin{displaymath}\langle Q^{\prime 2}_{xx}\rangle = \frac{1}{2\pi} \int_{-\inf...
P_{Q^{\prime}_{xx}} (\omega) {\rm d}\omega,
\end{displaymath} (18)

because $P_{Q_{xx}^{\prime}}(\omega)$ is an even function. The integral can be computed numerically because the power spectrum of the kinetic energy fluctuations decreases rapidly for $\omega > \omega_{0}$. A rough order-of-magnitude estimate can be found considering that $\omega_{0}$ is much larger than both $\omega_{\rm f}$ and $\lambda$:

 \begin{displaymath}\langle Q^{\prime 2}_{xx}\rangle \approx \frac{\langle {\cal T}^{\prime 2}\rangle }{\omega_{0}^{4}}\cdot
\end{displaymath} (19)

The mean square fluctuation of the quadrupole moment produced by the whole convection zone can be computed by co-adding the uncorrelated contributions of the single layers, as given by Eq. (18).

3 Fluctuations of the orbital elements in close binaries

In a close binary system, the fluctuations of the stellar quadrupole moment of the convective component produce fluctuating forces that act on the companion, which in turn induces a perturbation of the orbital motion. The variations of the orbital elements under the action of a random force have already been studied in other contexts and a general formalism has been developed in the framework of Hamiltonian mechanics (cf., e.g., Barge et al. 1982). For the present problem, however, we shall make use of the simplest formalism based on Lagrange planetary equations for the variation of the orbital elements in the formulation of Gauss (e.g., Roy 1978). In that approach the gradient of the perturbing function, the so-called perturbing force, associated with the fluctuating quadrupole moment is decomposed into three components: S, directed outward along the radius vector, T in the orbital plane at right angle to S, such that it makes an angle less than $90^{\circ}$ with the velocity vector, and W, perpendicular to the orbital plane and positive towards the north side of the plane. Starting from Eq. (1) and assuming a reference frame with the x-axis directed along the line joining the centers of the two stars from the convective component to its companion, the only component of the quadrupole moment tensor which contributes to the perturbing acceleration is Qxx and the perturbing force is:

\begin{displaymath}S = - \frac{9}{2} \, G \left( \frac{M + m}{M} \right) \frac{Q_{xx}}{r^{4}}

T = 0     (20)

W = 0

where M and m are the masses of the convective component and its companion star, respectively, and r is the orbital radius. For simplicity sake, we assume that the reference (i.e., unperturbed) orbit is circular with semi-major axis a.

Therefore, the only two Lagrange equations that do not reduce to an identity are:

 \begin{displaymath}\frac{{\rm d}\epsilon}{{\rm d}t} = 9\, G \left( \frac{M + m}{M} \right)
\frac{Q_{xx} (t)}{n a^{5}}
\end{displaymath} (21)


 \begin{displaymath}\frac{{\rm d}e}{{\rm d}t} = - \frac{9}{2} \,G \left( \frac{M + m}{M} \right)
\frac{Q_{xx}(t)}{na^{5}} \sin n t + t_{\rm v}
\end{displaymath} (22)

where $\epsilon$ is the mean longitude at the epoch, e is the eccentricity of the orbit, t is the time, $n \equiv 2\pi/P$ is the mean motion (P is the orbital period), and $t_{\rm v}$ is the tidal force which induces the circularization of the orbit (Zahn 1977). Following Zahn (1977), it is convenient to substitute Eq. (22) with its average taken over an orbital period:

 \begin{displaymath}\frac{{\rm d}e}{{\rm d}t} = - \frac{9}{2}\, G \left( \frac{M + m}{M} \right)
\frac{1}{na^{5}} f(t) -\frac{e}{\tau_{\rm c}}
\end{displaymath} (23)

where the average f(t) is defined as:

 \begin{displaymath}f(t) \equiv \frac{1}{P} \int_{t-P/2}^{t+P/2} Q_{xx}(t^{\prime}) \sin n t^{\prime} {\rm d} t^{\prime},
\end{displaymath} (24)

and the tidal circularization time $\tau_{\rm c}$ has been introduced. The solutions of the stochastic Eqs. (21) and (23) can be expressed in terms of the distribution functions $\psi_{\epsilon} (\epsilon, t)$ and $\psi_{\rm e} (e,t)$. The distribution function for $\epsilon$ is defined in such a way that the probability of finding $\epsilon$ in the interval $[\epsilon_{0}, \epsilon_{0}+{\rm d}\epsilon_{0}]$ at a given time t0 is: $\psi_{\epsilon} (\epsilon_{0}, t_{0}) {\rm d}\epsilon_{0}$; it also obeys a normalization condition at each time so that: $\int_{-\infty}^{+\infty} \psi_{\epsilon}
(\epsilon, t) {\rm d}\epsilon = 1$. Similar definition and normalization property are valid for the distribution function of the eccentricity $\psi_{\rm e} (e,t)$.

We shall consider long time intervals with respect to those characteristic of the fluctuations of Qxx, so that the stochastic variations of the perturbing force can be regarded as a Markovian process, i.e., with an autocorrelation function of the form: $\langle Q^{\prime}_{xx}(t_{1}) Q^{\prime}_{xx}(t_{2}) \rangle = \langle Q^{\prime \, 2}_{xx} \rangle \delta (t_{1} -t_{2})$, where $\delta$ is the Dirac delta function. By applying the theory outlined by, e.g., Rytov et al. (1987), we can write the Fokker-Planck equations ruling the probability density distributions $\psi_{\epsilon}$ and $\psi_{\rm e}$. We find:

 \begin{displaymath}\frac{\partial \psi_{\epsilon}}{\partial t} = \frac{1}{2} B
\frac{\partial^{2} \psi_{\epsilon}}{\partial \epsilon^{2}},
\end{displaymath} (25)

where the diffusion coefficient is:

 \begin{displaymath}B = 81 \, G^{2} \left( \frac{M+m}{M}\right)^{2} \frac{\langle Q^{\prime\, 2}_{xx} \rangle }{n^{2} a^{10}}
\end{displaymath} (26)

and the periodic boundary condition imposed on the probability density is: $\psi_{\epsilon} (\epsilon, t) = \psi_{\epsilon} (\epsilon + 2k\pi, t)$, with k a relative integer, and

 \begin{displaymath}\frac{\partial \psi_{\rm e}}{\partial t} =
...rac{1}{2}\, C
\frac{\partial^{2}\psi_{\rm e}}{\partial e^{2}},
\end{displaymath} (27)

where the diffusion coefficient is:

\begin{displaymath}C = \frac{81}{4}\, G^{2} \left( \frac{M+m}{M}\right)^{2} \frac{\langle f^{2} \rangle }{n^{2} a^{10}}\cdot
\end{displaymath} (28)

The boundary condition to be satisfied is $\psi_{\rm e} (e,t) =0$ for e<0, because the eccentricity can not be negative.

We are interested in the solution of Eq. (25) in a short time interval compared to the time needed by $\epsilon$ to change by $2\pi$by random walk fluctuations. Therefore, its solution is that typical of a diffusion process:

 \begin{displaymath}\psi_{\epsilon} (\epsilon, t) = \frac{1}{\sqrt{2\pi B t}} \exp \left(
-\frac{\epsilon^{2}}{2 B t}\right)\cdot
\end{displaymath} (29)

The diffusion coefficient B can be re-written introducing Kepler's third law into Eq. (26) and summing up the contributions of each layer of the convective envelope to the quadrupole moment fluctuations:

 \begin{displaymath}B = \frac{81 n^{2}}{M^{2}a^{4}} \sum_{L} \langle Q_{xx}^{\pri...
...{L} \frac{\langle {\cal T}^{\prime 2}\rangle }{\omega_{0}^{4}}
\end{displaymath} (30)

where L indicates the layers over which the summation is extended.

In the case of Eq. (27), we are interested in the stationary probability density for the residual eccentricity, that is attained after a sufficiently long time with respect to  $\tau_{\rm c}$. By assuming that the l.h.s. of Eq. (27) vanishes and imposing the boundary condition, we find:

 \begin{displaymath}\psi_{\rm e} (e, \infty) = \frac{1}{\sqrt{2\pi \langle e^{2} ...
...}} \exp \left( - \frac{e^{2}}{2 \langle e^{2} \rangle}\right),
\end{displaymath} (31)

where the mean square value of e is defined as:

 \begin{displaymath}\langle e^{2} \rangle = \frac{1}{2}\, C \, \tau_{\rm c}.
\end{displaymath} (32)

It is easier to obtain the mean square value of e from the power spectrum of its fluctuations $P_{\rm e} (\omega)$ derived from Eq. (23) (cf. Panter 1965):

\begin{displaymath}P_{\rm e} (\omega) = \frac{81 n^{2}}{4 M^{2}a^{4}}\frac{P_{\r...
...{\left( \frac{1}{\tau_{\rm c}^{2}} + \omega^{2}\right)} \cdot
\end{displaymath} (33)


 \begin{displaymath}\langle e^{2} \rangle = \frac{81 n^{2}}{4\pi M^{2} a^{4}} \in...
...c{1}{\tau_{\rm c}^{2}} + \omega^{2}\right)} \, {\rm d}\omega .
\end{displaymath} (34)

The denominator of Eq. (34) is very small for $\omega \ll \tau_{\rm c}^{-1}$and increases rapidly for $\omega \geq \tau_{\rm c}^{-1}$; since $\tau_{\rm c}^{-1}$is very small in comparison with the frequency of fluctuation of f, $P_{\rm f}(\omega)$ can be regarded as multiplied by a delta function centered at $\omega =0$. Therefore, the mean square value of e is given by:

 \begin{displaymath}\langle e^{2} \rangle = \frac{81 \tau_{\rm c} n^{2}}{8 M^{2} a^{4}}
P_{\rm f} (0).
\end{displaymath} (35)

In Appendix B it is shown that $P_{\rm f} (0) = (1/2) P_{Q^{\prime}_{xx}} (n)$; therefore, we arrive at the final result:

\begin{displaymath}\langle e^{2} \rangle = \frac{81 \tau_{\rm c} n^{2}}{16 M^{2} a^{4}} P_{Q^{\prime}_{xx}} (n).
\end{displaymath} (36)

The orbital frequency n is much larger than the characteristic frequency $\omega_{\rm f}$ in Eq. (16) (cf. Appendix A), hence we can neglect it in the expression of $P_{Q^{\prime}_{xx}}$ and obtain:

 \begin{displaymath}\langle e^{2} \rangle = \frac{18 \tau_{\rm c}}{M^{2} a^{4}}
\frac{P_{\cal T}(n)}{(n^{2} + \lambda^{2})} \cdot
\end{displaymath} (37)

It is important to note that Eq. (37) gives the contribution to the residual eccentricity produced by a single convective layer. The residual eccentricity produced by the entire convective envelope can be computed by summing up the power spectra of the kinetic energy of the single layers in the hypothesis of uncorrelated fluctuations:

 \begin{displaymath}\langle e^{2} \rangle = \frac{18 \tau_{\rm c}}{ M^{2} a^{4} (n^{2} + \lambda^{2})}
{\sum_{L} P_{\cal T}(n)},
\end{displaymath} (38)

where again L denotes the layers over which the summation is extended.

4 Application to binaries containing millisecond pulsars

In this section, an application of the theory outlined above will be presented dealing with the residual eccentricity of binaries containing millisecond pulsars. Moreover, we shall compute the order of magnitude of the fluctuations of the mean longitude at the epoch and show that they are too small to be observed.

4.1 The residual eccentricity

As reviewed by Verbunt (1990, 1993) and van den Heuvel (1992), binary millisecond pulsars are believed to have been spun-up by mass transfer and their progenitors should have been X-ray binaries. According to Phinney (1992) and Phinney & Kulkarni (1994), the observed small eccentricities of their orbits are fossil residuals of their spin-up phase.

The evolutionary channels that may produce binary millisecond pulsars have been recently reviewed by Taam et al. (2000) and we refer the reader to that paper for more details. For our purposes it is interesting to consider the systems that may be formed through two of the proposed channels: a) late evolution of a low-mass secondary, in which the donor star during the mass transfer phase was a red giant with a degenerate helium core and ended its evolution as a helium white dwarf (initial mass $M< 2.25~M_{\odot}$); b) early massive Case B evolution, in which the donor star was massive enough to ignite helium under non-degenerate conditions; if the initial mass ratio of the binary was smaller than a critical value (which depends on the initial orbital separation, see King & Ritter 1999), the final phase of mass transfer took place after the end of the central helium burning, when the star had approached the AGB branch. At this stage it has an extended outer convective envelope. The donor star ends its evolution as a CO white dwarf, after having depleted its envelope (cf. Tauris et al. 2000).

In both scenarios the final phase of mass transfer is characterized by a secondary with a deep convective envelope. The matter flowing through the inner Lagrangian point to the primary pulsar is responsible for its spin-up and X-ray emission so that the system is observable as a low-mass X-ray binary. That phase comes to an end when the mass of the secondary's envelope falls below ${\approx} 0.01~M_{\odot}$ and the star starts to shrink and detaches from its Roche lobe. It is during this stage of the secondary evolution that the residual eccentricity is produced, according to the model proposed by Phinney (1992). The shrinking secondary still retains a convective envelope which is responsible for the fluctuations of the gravitational quadrupole moment of the star, thus inducing a residual eccentricity of the orbit. The contraction occurs on a time scale ranging from ${\sim} 7\times 10^{5}$ to ${\sim} 10^{7}$ yr in the case of a degenerate helium core secondary, thus, when the secondary radius is close to the radius of its Roche lobe $R_{\rm L}$, the circularization time $\tau_{\rm c} \sim 10^{4}$ yr is much shorter than the evolutionary time and the mean square value of the eccentricity is given by Eq. (38). However, once the radius of the secondary R becomes ${\sim} 0.5\,R_{\rm L}$, the timescale for pumping and dumping the eccentricity becomes longer than the time scale on which the envelope is contracting, because $\tau_{\rm c} \propto R^{-8}$. The eccentricity is then "frozen" at its value when $R \approx 0.5\, R_{\rm L}$.

The radius of the Roche lobe is given with a good approximation by: $R_{\rm L}= 0.49 a q^{2/3}/[0.6 q^{2/3} + \ln (1 + q^{1/3})]$, where a is the semi-major axis of the orbit and q=M/m is the mass ratio (Eggleton 1983). For simplicity, we shall assume that the mass of the pulsar is $m=1.4~M_{\odot}$in all systems. The radius of the secondary star is primarily a function of its core mass (cf. Refsdal & Weigert 1970; Ritter 1999) and at the end of the mass transfer phase it equals the maximum radius corresponding to the mass of its core, that is the white dwarf companion of the millisecond pulsar. Refsdal & Weigert (1971) gave a simple relationship between the maximum radius of the secondary and its core mass: $R_{\rm max}/R_{\odot} = 10^{4.1} (M/M_{\odot})^{5.3}$, valid for stars of solar chemical composition and neglecting the small difference between the mass of the secondary and that of its helium core ($\leq$5%). Since the maximum radius is equal to the Roche lobe radius at the end of mass transfer, we can express the semi-major axis and the orbital period as a function of M. More accurate expressions of the mass-period relation are given by King & Ritter (1999) and Tauris & Savonije (1999).

Such a mass-period relation is a consequence of the proposed evolutionary scenario and its validity may be used to distinguish between binary pulsars that were formed through stable mass transfer from a donor star with a core-envelope structure, and binary pulsars that had a different origin (cf. Phinney & Kulkarni 1994; King & Ritter 1999; Taam et al. 2000).

We have applied this criterium to the systems listed by Taam et al. (2000) considering only those systems with a measured residual eccentricity or a significant upper limit for it (which excludes J1435-60, J1232-6501, J1453-58 and J1810-2005). The orbital period P of these 30 systems is plotted versus the most probable value of their secondary's mass M in Fig. 1, together with the theoretical P-M relationships derived by King & Ritter (1999) and Refsdal & Weigert (1971), respectively.

\par\includegraphics[width=8cm,clip]{h2694f1.ps}\end{figure} Figure 1: The orbital period P of the millisecond binary pulsars listed in Table 1 of Taam et al. (2000) versus their secondary's mass, together with the theoretical P-M relationships given by King & Ritter (1999) ( solid line) and Refsdal & Weigert (1971) ( dashed line). The filled circles with error bars give the most probable values of their mass (corresponding to the median inclination $i=60^{\circ }$), while the horizontal error bars indicate the 90% confidence intervals of their masses. Open diamonds indicate systems that do not agree with the theoretical P-M relationship or that have a residual eccentricity larger than 0.01, which is not compatible with the evolutionary scenario discussed in the text.
Open with DEXTER

It appears that B1957+20, J2051-0827, B1831-00, J1603-7202, J2145-0750, J0621+1002, J1022+1001 and B0655+64 (open symbols) are not consistent with the proposed evolutionary scenario. Moreover, we disregard also J1904+04 whose orbital eccentricity is 0.04, too large to be explained by the mechanism proposed by Phinney (1992).

It is important to note that systems with an orbital period smaller than about 1-2 days, although in several cases appear to agree with the theoretical P-M relationship, should have had a different evolution because the mass transfer started when their secondaries were still on the main sequence (Ergma et al. 1998; Taam et al. 2000).

In conclusion, the comparison with the theoretical P-M relationship and the requirement that P > 2 days lead us to exclude 13 systems. For the remaining 17 systems we can compare, as done by Phinney & Kulkarni (1994), the observed residual eccentricity with the mean square value predicted according to the evolutionary scenario outlined above. In the present paper we do not perform a detailed comparison, but make an illustrative application in order to show how the approach developed in Sects. 2 and 3 is applicable to a stellar model. We adopt the models calculated by Driebe et al. (1998) for the evolution of secondaries with a degenerate He core from the red giant branch to the white dwarf sequence and final mass ranging from $0.179~M_{\odot}$ to 0.414 $~M_{\odot}$. The luminosity L and the radius Rare derived from their evolutionary tracks, whereas the mass of the convective envelope $M_{\rm env}$ and the radius at the base of the convective envelope $r_{\rm b}$ are assumed to be constant for a given model. The values of $M_{\rm env}$ and $r_{\rm b}$ are listed in Table 1.

Table 1: Model parameters for the secondary components.
M $M_{\rm env}$ $r_{\rm b}$
($M_{\odot}$) ( $10^{-3}~M_{\odot}$) ($R_{\odot}$)
0.179 48.289 0.0261
0.195 30.768 0.0246
0.234 13.098 0.0219
0.259 7.232 0.0208
0.300 4.746 0.0192
0.331 3.744 0.0181
0.414 2.175 0.0159

Our model of the convective envelope of the donor star assumes that it is adiabatic and has a negligible mass in comparison with the core mass. For a perfect gas with a ratio of the specific heats $\gamma =5/3$, the density $\rho$ is given by:

 \begin{displaymath}\rho(r) = \rho_{0} \left(\frac{R-r}{r}\right)^{3/2}
\end{displaymath} (39)


\begin{displaymath}\rho_{0}= \frac{M_{\rm env}}{4\pi R^{3} \int_{z_{\rm b}}^{1} z^{1/2} (1-z)^{3/2} {\rm d}z},
\end{displaymath} (40)

r is the distance from the centre of the star, and $z_{\rm b} \equiv
r_{\rm b}/R$.

The pressure scale height is given by:

 \begin{displaymath}H_{\rm p} = \frac{2 (R -r) r}{5R}
\end{displaymath} (41)

and the mixing length $\ell$ is assumed to be proportional to $H_{\rm p}$: $\ell = \alpha H_{\rm p}$. The convective velocity is computed from the mixing-length theory by assuming that the convective flux is given by $F_{\rm c} = (10/\alpha) \rho u^{3}$and that the energy transport through the envelope is due only to convective motions (cf., e.g., Kippenhahn & Weigert 1990). This yields:

\begin{displaymath}u = \left(\frac{\alpha L}{40\pi r^{2} \rho}\right)^{1/3}.
\end{displaymath} (42)

Therefore, the mean square value of the fluctuation of the kinetic energy in a layer of thickness $\ell$ is:

\begin{displaymath}\langle{\cal T}^{\prime 2} \rangle \simeq r^{2} \overline{\rho^{2}} u^{4} \ell^{4},
\end{displaymath} (43)

where $\overline{\rho^{2}}$ is the mean square density in the layer and $\omega_{0} = u/\ell$ the turn-over frequency. The diffusion time for the envelope is $\tau_{\rm d}=R^{2}
/(\ell_{\rm b} u_{\rm b})$, where $\ell_{\rm b}$ is the mixing length and $u_{\rm b}$ the convective velocity at the base of the envelope $r=r_{\rm b}$. The power spectrum of the kinetic energy fluctuations can be evaluated from Eq. (17) for any given layer. Then the power spectrum of the fluctuations for the whole convective envelope is computed by summing up the contributions of each layer, the radial extension of which is assumed to be equal to its mixing length. According to Verbunt & Phinney (1995), the circularization rate is given by:
$\displaystyle \tau_{\rm c}^{-1}$=$\displaystyle g \left(\frac{L}{M_{\rm env} R^{2}}\right)^{1/3}\left(\frac{M_{\rm env}}{M}\right)$
$\displaystyle \times\left(\frac{m+M}{M}\right) \left(\frac{m}{M}\right)
\left(\frac{R}{a}\right)^{8}$ (44)

where g is a dimensionless parameter coming from the theory of stellar tides, m is the mass of the pulsar and the other symbols have already been defined.

The final residual eccentricity of the binary comes out as the result of random excitations and dumpings due to the fluctuating gravitational field of the secondary and to tidal dissipation, respectively. The contribution to the final eccentricity produced during a small time interval $[t, \; t+\Delta t]$ (with $\Delta t < \tau_{\rm c}$) is:

\begin{displaymath}\Delta \langle e^{2}_{\rm res} \rangle \simeq
\langle e^{2} \rangle (t) \exp[-2\theta(t)] \Delta t/\tau_{\rm c}(t),
\end{displaymath} (45)

where $\langle e^{2} \rangle (t)$ is the value of the mean square eccentricity computed by Eq. (38) with the instantaneous values of $\tau_{\rm c}$ and $P_{\cal T}$, and

 \begin{displaymath}\theta(t) = \int_{t}^{t_{\rm s}} \frac{{\rm d}t^{\prime}}{\tau_{\rm c}(t^{\prime})},
\end{displaymath} (46)

where $t_{\rm s}$ is the age of the system, is the variation of the logarithm of the eccentricity in the time interval $[t, t_{\rm s}]$ (cf. Verbunt & Phinney 1995). We assume $t_{\rm s}=10$ Gyr for $M=0.179~M_{\odot}$, and $t_{\rm s}=7$ Gyr in all of the other cases.

Therefore, the mean square value of the residual eccentricity is approximately given by:

 \begin{displaymath}\langle e_{\rm res}^{2} \rangle \simeq \int_{0}^{t_{\rm s}} \...
...{2} \rangle (t)}{\tau_{\rm c}(t)} \exp [-2\theta(t)] {\rm d}t.
\end{displaymath} (47)

The final residual eccentricity has been computed for the seven models published by Driebe et al. (1998) by deriving the orbital period and the semi-major axis of the orbit from the mass-period relation of King & Ritter (1999) and assuming that the mass of the pulsar is $m=1.4~M_{\odot}$. The mixing length parameter $\alpha$ was assumed to be 1.7 for consistency with the calculations of Driebe et al. (1998), and the tidal parameter g=0.8, in agreement with the prescription of the tidal theory (cf. Verbunt & Phinney 1995). In order to compute the integrals in Eqs. (46) and (47), we took into account also the phases of rapid radius and luminosity variations associated with the thermal instabilities of the hydrogen burning energy source for $0.21 \protect\raisebox{-0.5ex}{$\:\stackrel{\textstyle <}
{\sim}\:$ }M/M_{\odot} \protect\raisebox{-0.5ex}{$\:\stackrel{\textstyle <}
{\sim}\:$ }0.3 $ (cf. Driebe et al. 1999). Such phases turn out to be important because the radius varies on a time scale comparable with $\tau_{\rm c}$, thus giving rise to "impulsive'' contributions to the residual eccentricity that can be of the order of ${\sim} 10$% of its value. The theoretically predicted $\langle e_{\rm res}^{2} \rangle$ are compared with the observed eccentricities in Fig. 2, where a quadratic fit to the theoretical results is also displayed.

4.2 Fluctuations of the mean longitude at the epoch

The coefficient B appearing in the distribution of the mean longitude at the epoch (Eq. (29)) can be computed according to Eq. (30) using Eqs. (16) and (18), where the contributions coming from the different convective layers of the star are to be added. The contribution of each layer can be evaluated numerically for the secondaries considered above and we find that B is in any case extremely small ( ${\sim} 3\times 10^{-109} < B < \,\sim 3 \times 10^{-121}$ yr-1 for 2 < P < 1000 days). Therefore, there is no possibility of observing the fluctuations of the mean longitude $\epsilon$ during the phase when the system was a LMXB since the accuracy with which we can measure the variations of $\epsilon$ is of the order of 10-4-10-3 rad for eclipsing systems observed over several tens of years.

\par\includegraphics[width=7.8cm,clip]{h2694f2.ps}\end{figure} Figure 2: The observed eccentricities $e_{\rm res}$versus the orbital period P for our sample of binary millisecond pulsars (filled circles) and the theoretical values of $\langle e_{\rm res}^{2} \rangle^{1/2}$ calculated for the seven models of Driebe et al. (1999) (open circles). The solid line is a quadratic least square fit to the theoretical values.
Open with DEXTER

5 Discussion and conclusion

The application of the virial theorem allows us to derive an expression which relates the fluctuations of the gravitational quadrupole moment of the star with the kinetic energy fluctuations of the convective motions. This provides a rigorous theoretical base to the dissipation-fluctuation theorem proposed by Phinney (1992) and supports the interpretation he proposed for the residual eccentricities of binaries containing millisecond pulsars. Our illustrative calculations show that the mean square value of the residual eccentricity is a function of the secondary's mass only, i.e., of the orbital period, given the mass-period relationship for binary millisecond pulsars. The agreement between our theoretical predictions and the observations is good for orbital periods in the range 2-10 d, but significant deviations of a factor of ${\sim} 10$ appear for systems with periods longer than 20 days. This discrepancy may be possibly ascribed to the limitations of the adopted secondary star models. Specifically, Driebe's et al. models are characterized by maximum radii which are significantly smaller than the Roche lobe radii of long period systems. This may be due to the assumptions made on the mass loss rate during the evolution of the secondaries, which was not derived from a detailed physical model. As a matter of fact, the evolutionary calculations by Tauris & Savonije (1999), who included a more sophisticated treatment of mass loss, yield larger maximum radii for the secondaries. It is interesting to note that our assumption of constant envelope mass is not critical because the dependencies on $M_{\rm env}$ cancels each other in the expression for $\langle e^{2} \rangle$. Such a conclusion is also confirmed by some exploratory computation with ${\dot{M}}_{\rm env}\propto L$. Moreover, the adopted values for $r_{\rm b}$ do not significantly influence our results, provided that $R\gg r_{\rm b}$ in the relevant evolutionary phases.

Assuming that the tidal circularization process is properly described, the distribution of the residual eccentricities may be used to test models of the secondary's evolution and to probe the properties of stellar convection, in particular the kinetic energy spectrum, in the phases following the end of mass transfer. The approach presented above is simply illustrative, also in view of the limitations of the present evolutionary models. The development of detailed models of the secondary's evolution, which are capable of following its radius and luminosity variations in a realistic way during all of the relevant phases, will allow us to make accurate predictions by numerically solving the initial value problem posed by Eq. (31) with a time varying coefficient C and computing the change of the distribution $\psi_{\rm e} (e,t)$ versus time. We wish to apply such a more rigorous approach in a forthcoming paper.

The proposed method to compute quadrupole moment fluctuations is interesting also for other problems that involve binary systems with a convective component and a circularization time $\tau_{\rm c}$ long enough to allow a random walk growth of the residual eccentricity to significant values. The recently discovered extrasolar planets may be relevant in this respect. In those systems, the much smaller mass of the planet, with respect to the star, makes  $\tau_{\rm c}$ very long, thus allowing the residual eccentricity to grow. However, an exploratory calculation we performed by assuming the central star to have the same structure of the present Sun and a 4.0 d orbital period for the planet, indicates that the attained residual eccentricities are in the range $10^{-8}{-}5\times 10^{-7}$, even if the tidal circularization time may well exceed the star main sequence life-time (cf. Trilling 2000). However, the situation may be significantly different if the planet is orbiting a red giant because in that case the amplitude of the convective fluctuations is much more pronounced and effective.

The authors wish to thank an anonymous Referee for a careful reading of the manuscript and stimulating comments. Research on stellar physics at Catania Astrophysical Observatory and at the Dept. of Physics and Astronomy of Catania University is funded by MURST (Ministero dell'Università e della Ricerca Scientifica e Tecnologica) and Regione Sicilia, whose financial support is gratefully acknowledged.

Appendix A

In this Appendix the relation between the variation of the tensor $\Xi_{ik}$ and the gravitational quadrupole moment Qik is discussed. The variation of the gravitational energy tensor due to the density fluctuations can be written as:

$\displaystyle {\cal G}_{ik}^{\prime} = -\frac{G}{2}$  
$\displaystyle \times \int_{V}\int_{V}
\left[ \rho_{0}({\vec x}) \rho^{\prime}({...
...}({\vec x},t)
+ \rho^{\prime}({\vec x})\rho^{\prime}({\vec x}^{\prime}) \right]$  
$\displaystyle \times \frac{(x_{i}-x_{i}^{\prime})(x_{k} -x_{k}^{\prime})}{\vert{\vec x} - {\vec x^{\prime} }\vert^{3}} {\rm d}{\vec x} {\rm d}{\vec x^{\prime}}.$ (48)

The terms containing the products $\rho_{0}({\vec x})\rho^{\prime}({\vec x}^{\prime},t)$ and $\rho_{0}({\vec x}^{\prime})\rho^{\prime}({\vec x},t)$are averaged to zero by the volume integration and only the term containing $\rho^{\prime}({\vec x},t)\rho^{\prime}({\vec x}^{\prime},t)$ will survive. Applying the mean value theorem, it can be written as:
$\displaystyle {\cal G}_{ik}^{\prime} = -\frac{G}{2L^{3}_{ik}} \bigg\langle \left( \frac{\rho^{\prime}}{\rho_{0}} \right)^{2} \bigg\rangle$  
$\displaystyle \times \int_{V} \int_{V} \rho({\vec x},t)
\rho({\vec x}^{\prime},...
...{i}^{\prime}) (x_{k} - x_{k}^{\prime})
{\rm d}{\vec x} {\rm d}{\vec x}^{\prime}$ (49)

where $\bigg\langle \left( \frac{\rho^{\prime}}{\rho_{0}} \right)^{2} \bigg\rangle$ indicates the mean squared density fluctuation averaged over the convection zone and the average length Lik depends on the distribution of mass inside the stellar configuration and the spatial direction indexes i and k. If the origin of the reference frame is fixed at the barycentre of the star, it yields:

\begin{displaymath}{\cal G}^{\prime}_{ik} = -\frac{GM}{L^{3}_{ik}} \bigg\langle ...
...{\rho^{\prime}}{\rho_{0}} \right)^{2} \bigg\rangle I_{ik}(t),
\end{displaymath} (50)

where the time dependence of the inertia tensor, coming from Eq. (49), has been explicitly noted.

We define the characteristic frequency $\omega_{\rm f}$ as:

\begin{displaymath}\omega_{\rm f}^{2} \equiv \frac{GM}{L^{3}_{ik}} \bigg\langle \left( \frac{\rho^{\prime}}{\rho_{0}} \right)^{2} \bigg\rangle,
\end{displaymath} (51)

assuming that it is independent of the time and the spatial directions for a stationary and spherically symmetric fluctuating density field. Therefore, we can write:

\begin{displaymath}{\cal G}^{\prime}_{ik} = -\omega_{\rm f}^{2} I_{ik}(t).
\end{displaymath} (52)

The time variation of $\Xi_{ik}$ then immediately follows:

\begin{displaymath}\Xi_{ik}^{\prime} = - \omega_{\rm f}^{2} Q_{ik}^{\prime}.
\end{displaymath} (53)

It is interesting to estimate the order of magnitude of the characteristic frequency $\omega_{\rm f}$ for a giant burning hydrogen in a thin shell at the base of its convective envelope, which may be relevant to the case of the progenitors of binary millisecond pulsars, and for a main sequence star like the Sun. In the first case we can assume: $M = 0.35~M_{\odot}$, $R= 50~R_{\odot}$ and $\bigg\langle \left( \frac{\rho^{\prime}}{\rho_{0}} \right)^{2} \bigg\rangle \approx \left( \frac{u}{c_{\rm s}}\right)^{4} \sim 10^{-6}$, where u is the average convective velocity and $c_{\rm s}$ the average sound speed; this yields: $\omega^{2}_{\rm f} = 1.1 \times 10^{-18}$ s-2. In the case of the present Sun, $\bigg\langle \left( \frac{\rho^{\prime}}{\rho_{0}} \right)^{2} \bigg\rangle \sim 10^{-8}$, which yields: $\omega^{2}_{\rm f} = 3.9\times 10^{-15}$ s-2.

Appendix B

In this Appendix we shall compute the Fourier transform of the function f(t) defined in Eq. (24). Introducing Euler identity, we can write:

\begin{displaymath}f(t)= \frac{1}{P} \int_{t-P/2}^{t+P/2} Q_{xx}^{\prime}(t^{\pr...
...\prime})} - \exp {(-i n t^{\prime})}}{2 i} {\rm d} t^{\prime}
\end{displaymath} (54)

and, introducing the Fourier transform of $Q_{xx}^{\prime}(t)$:
f(t)=$\displaystyle \frac{1}{P} \int_{t-P/2}^{t+P/2} \frac{1}{2\pi} \int_{-\infty}^{\infty}
\tilde{Q}_{xx}^{\prime} (\omega)
\exp {i\omega t^{\prime}}$
$\displaystyle \times \frac{\exp {(i n t^{\prime})} - \exp {(-i n t^{\prime})}}{2 i}
{\rm d} t^{\prime} {\rm d}\omega.$ (55)

Collecting together the arguments of the exponential functions and performing the integration over the variable $t^{\prime}$, we find:
f(t)=$\displaystyle \frac{1}{2\pi} \int_{-\infty}^{\infty} \frac{1}{P} \left[\tilde{Q}_{xx}^{\prime}(\omega-n)
-\tilde{Q}_{xx}^{\prime}(\omega + n) \right]$
$\displaystyle \times \frac{\sin \frac{\omega P}{2}}{i\omega} \exp i \omega t \, {\rm d}\omega .$ (56)

Therefore, the Fourier transform of f is:

\begin{displaymath}\tilde{f}(\omega) = \frac{1}{P} \left[\tilde{Q}_{xx}^{\prime}...
...e}(\omega + n) \right] \frac{\sin \frac{\omega P}{2}}{i\omega}
\end{displaymath} (57)

and, since $\tilde{Q}_{xx}(\omega) = \tilde{Q}_{xx}^{*}
(-\omega)$ and $\Re \tilde{Q}_{xx}(\omega) = \Im \tilde{Q}_{xx}(\omega)$ in a statistical sense (with the asterisk denoting complex conjugation and $\Re$ and $\Im$ the real and imaginary part, respectively), the relationship $P_{\rm f} (0) = (1/2) P_{Q^{\prime}_{xx}} (n)$ follows.



Copyright ESO 2001