text:
# The evolution of naked helium stars with a neutron-star companion in close binary systems

## 1 Introduction

A binary pulsar with a neutron-star or a heavy white-dwarf companion has long been considered to originate from a helium star in a binary system with a neutron-star companion. The latter system is a descendant of a high-mass X-ray binary (HMXB), in which the companion of the neutron star loses its mass through wind mass loss or a mass-transfer phase, exposing its helium core (e.g. Bhattacharya & van den Heuvel 1991). Although the existence of such a system is only found in Cyg X-3 (van Kerkwijk et al. 1992, 1996), a detailed study of the evolution of helium stars in binary systems with a compact companion is important as the systems form the bridge between the evolution of X-ray binaries and the formation of double compact-object binaries.

The evolution of a helium star in a binary system has been studied e.g. by Savonije, de Kool & van den Heuvel (1986) who evolved a 0.6 $`\mathrm{M}_{\mathrm{}}`$ non-degenerate helium star with a 1.3 $`\mathrm{M}_{\mathrm{}}`$ compact companion. With a short orbital period of $`P=37^\mathrm{m}`$, Roche-lobe overflow (RLOF) takes place during helium core burning. Ergma & Fedorova (1990) evolved helium stars with masses of 0.5, 0.766, and 1 $`\mathrm{M}_{\mathrm{}}`$. As the companion star, they took white dwarfs with the same range of masses in a combination such that the systems have a mass ratio of $`0.5M_{\mathrm{He}}/M_{\mathrm{WD}}2`$. The periods are also so short, of $`26\stackrel{m}{.}2`$$`62\stackrel{m}{.}6`$, that the helium stars transfer mass to the companion while they are still burning helium in the centre. A similar work was also carried out by Tutukov & Fedorova (1990). A study of systems with more massive helium stars in wider orbits has been carried out by e.g. Delgado & Thomas (1981) who considered helium stars with masses of 2, 2.7, 3.3, and 4 $`\mathrm{M}_{\mathrm{}}`$ with a massive main-sequence companion. The helium star fills its Roche lobe after helium is exhausted in the core. Habets (1986a) evolved a 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium star with a 17 $`\mathrm{M}_{\mathrm{}}`$ main-sequence companion in a wide orbit ($`P=20\stackrel{d}{.}25`$) such that mass-transfer phase occurs during first carbon burning convective shell.

In this work we study the evolution of helium stars with masses of 1.5 – 6.7 $`\mathrm{M}_{\mathrm{}}`$ with a 1.4 $`\mathrm{M}_{\mathrm{}}`$ neutron-star companion. Such systems are assumed to be the remnants of Be/X-ray binaries with masses in the range of 8 – 20 $`\mathrm{M}_{\mathrm{}}`$ which underwent mass transfer as case B (RLOF is initiated during hydrogen shell burning) or case C (during helium shell burning). As the result of the large mass ratio, mass transfer from the Be star to the neutron star is dynamically unstable and the two components are embedded in a common envelope (CE), leading to the spiraling-in of the neutron star in the envelope of the Be star. We assume that the system survived the common envelope and spiral-in process. We are interested to investigate the final fate of the systems (whether they will become white-dwarf/neutron-star or double neutron-star binaries), the type of supernovae (SN) they might produce, as well as to study which systems are stable to RLOF. A study on the similar range of mass (2 – 6 $`\mathrm{M}_{\mathrm{}}`$ helium stars) has been carried out by Avila Reese (1993). With a Roche radius of 0.6 and 0.7 $`\mathrm{R}_{\mathrm{}}`$, in that study RLOF takes place after helium core burning is terminated and the calculations were done up to the ignition of carbon in the centre. Apart from the slightly wider range of mass in our study compared to the latter work and a larger range of orbital periods used in our work, we also are able to follow the evolution of the helium stars to more advanced evolutionary stage, i.e. beyond the carbon ignition, which enables more detailed conclusions on the fate of the systems. We take 6.7 $`\mathrm{M}_{\mathrm{}}`$ as the upper limit of our calculation because more massive helium stars undergo a dynamically unstable mass transfer.

In Sect. 2 we describe the possibility for the formation channel of a helium star in a binary system with a neutron-star companion, and some basic assumptions used in the calculation of the orbital evolution. The computational code, the input parameters used in the code, and the calculation that constrains the initial mass of the helium star and the initial period are described in Sect. 3. We discuss the results in two sections; Sect. 4 for the cases in which RLOF takes place during helium core burning and Sect. 5 for that during helium shell burning. Our conclusions are given in Sect. 6.

## 2 The evolution of helium star in binary system

### 2.1 Formation of a helium star with a neutron-star companion in a binary system

An HMXB consists of a massive OB-type star and a compact object as the X-ray source. There are two classes of HMXB, i.e.: (i) OB-supergiant (standard massive) X-ray binaries, in which the X-ray source is powered by wind accretion from its almost Roche-lobe filling OB-supergiant companion, and (ii) Be/X-ray binaries, in which the X-ray source is powered by accretion of material when the compact object moves through the dense disc of its unevolved Be-star companion at periastron passage. The first class is found in close, almost circular orbits with companion stars more massive than 20 $`\mathrm{M}_{\mathrm{}}`$; while the latter in wide, eccentric orbits, with 8 – 20 $`\mathrm{M}_{\mathrm{}}`$ companion stars.

As the companion starts to fill its Roche lobe, a runaway mass transfer to the compact object will follow due to the high mass ratio of the system, leading to a CE and spiral-in phase. The outcome of this phase is a system consisting of the helium core of the companion and a neutron star. The helium star evolves and, in turn, might fill its Roche lobe. If the helium star is more massive than 2.2 $`\mathrm{M}_{\mathrm{}}`$ (Habets 1986a) and it does not further lose mass, it will collapse and produce a neutron star; and the system becomes a double neutron-star binary. We refer to Bhattacharya & van den Heuvel (1991) for review on the so-called ”standard scenario”111An evolutionary channel for the formation of double neutron-star binary similar to the standard scenario was also proposed by Tutukov & Yungelson (1973). They predicted a helium-star mass of 2 $`\mathrm{M}_{\mathrm{}}`$ as the limiting mass for white-dwarf/neutron-star formation..

An alternative to the standard model was proposed by Brown (1995). This new evolutionary scenario suggests that a double neutron-star binary is a descendant of a double helium-star binary, which is produced from a binary system with components of very similar masses. This twin-mass model was proposed because of the argument that a neutron star would be able to accrete matter at high rate during the CE and spiral-in process and will then become a black hole (e.g. Chevalier 1993), which would exclude the possibility that double neutron stars would result from the CE evolution of HMXB (the standard scenario). Nevertheless, whether or not a neutron star in a CE will indeed go to highly super-Eddington accretion is still very uncertain (see e.g. Chevalier 1996). Therefore, there is still room for the formation of systems consisting of a helium star with a neutron-star companion like in the standard scenario. Strong accretion onto the neutron star can be avoided, for example, if the donor star is very extended (Chevalier 1993) and if rotational effects are taken into account (Chevalier 1996).

### 2.2 Orbital evolution

The orbital angular momentum of a binary system which contains a helium star of mass $`M_{\mathrm{He}}`$ and a neutron star of mass $`M_{\mathrm{NS}}`$ is given by

$`J_{\mathrm{orb}}^2`$ $`=`$ $`{\displaystyle \frac{M_{\mathrm{He}}^2M_{\mathrm{NS}}^2}{M}}aG`$ (1)
where $`M=M_{\mathrm{He}}+M_{\mathrm{NS}}`$ is the total mass of the system, $`a`$ is the orbital separation, and $`G`$ is the constant of gravity. The rate of change in orbital separation of the system is expressed by

$`{\displaystyle \frac{\dot{a}}{a}}`$ $`=`$ $`2{\displaystyle \frac{\dot{J}_{\mathrm{orb}}}{J_{\mathrm{orb}}}}2\left({\displaystyle \frac{\dot{M}_{\mathrm{He}}}{M_{\mathrm{He}}}}+{\displaystyle \frac{\dot{M}_{\mathrm{NS}}}{M_{\mathrm{NS}}}}\right)+{\displaystyle \frac{\dot{M}}{M}}`$ (2)
The total change in orbital angular momentum is assumed to be affected by gravitational-wave radiation and by the loss of mass with angular momentum from the system, i.e.

$`{\displaystyle \frac{\dot{J}_{\mathrm{orb}}}{J_{\mathrm{orb}}}}`$ $`=`$ $`{\displaystyle \frac{\dot{J}_{\mathrm{gwr}}}{J_{\mathrm{orb}}}}+{\displaystyle \frac{\dot{J}_{\mathrm{ml}}}{J_{\mathrm{orb}}}}`$ (3)
where the rate of change of orbital angular momentum due to gravitational-wave radiation is given by (Landau & Lifshitz 1958)

$`{\displaystyle \frac{\dot{J}_{\mathrm{gwr}}}{J_{\mathrm{orb}}}}`$ $`=`$ $`{\displaystyle \frac{32G^3}{5c^5}}{\displaystyle \frac{M_{\mathrm{He}}M_{\mathrm{NS}}M}{a^4}}`$ (4)
with $`c`$ the speed of light in vacuum.

The system is assumed to evolve non-conservatively. Mass loss from the system causes loss of orbital angular momentum which is taken into account as (van den Heuvel 1994; Soberman, Phinney & van den Heuvel 1997)

$`{\displaystyle \frac{\dot{J}_{\mathrm{ml}}}{J_{\mathrm{orb}}}}`$ $`=`$ $`{\displaystyle \frac{\alpha +\beta q^2}{1+q}}{\displaystyle \frac{\dot{M}_{\mathrm{He}}}{M_{\mathrm{He}}}}`$ (5)
where $`q=M_{\mathrm{He}}/M_{\mathrm{NS}}`$ is the mass ratio; $`\alpha `$ is the fraction of mass lost from the helium star in the form of fast isotropic wind, and $`\beta `$ is the fraction of mass ejected isotropically from the vicinity of the neutron star. During the detached phase, we took $`\alpha =1`$ and it is assumed that the neutron star does not accrete matter from the stellar wind ($`\beta =0`$); all matter released from the system in this case carries the specific orbital angular momentum of the helium star. The stellar wind mass loss is approximated by eq. (2) in Wellstein & Langer (1999, cf. Hamann, Schönberner & Heber 1982; Hamann, Koesterke & Wessolowski 1995) multiplied by a factor of 0.5 to provide a better fit with the most recent observed mass-loss rate of Wolf-Rayet stars (Nugis & Lamers 2000; Nelemans & van den Heuvel 2001), i.e.

$`\dot{M}_{\mathrm{He},\mathrm{wind}}`$ $`=`$ $`\{\begin{array}{cc}\mathrm{2.8\hspace{0.17em}10}^{13}(L/\mathrm{L}_{\mathrm{}})^{1.5},\hfill & \mathrm{log}L/\mathrm{L}_{\mathrm{}}4.5\hfill \\ \mathrm{4.0\hspace{0.17em}10}^{37}(L/\mathrm{L}_{\mathrm{}})^{6.8},\hfill & \mathrm{log}L/\mathrm{L}_{\mathrm{}}<4.5\hfill \end{array}`$ (8)
This means our adopted mass loss rate in the upper luminosity range is one quarter of the Hamann et al. (1995) rate. During RLOF, the helium star transfers mass with a mass-transfer rate which is defined as

$`\dot{M}_{\mathrm{He},\mathrm{RLOF}}`$ $`=`$ $`10^3\times \mathrm{max}[0,\left(\mathrm{ln}{\displaystyle \frac{R}{R_\mathrm{L}}}\right)^3]`$ (9)
where $`R_\mathrm{L}`$ is the helium star’s Roche radius given by Eggleton (1983) as

$`R_\mathrm{L}`$ $`=`$ $`{\displaystyle \frac{0.49a}{0.6+q^{2/3}\mathrm{ln}(1+q^{1/3})}}`$ (10)
During the mass-transfer phase (even on a nuclear timescale) the wind mass loss contributes at most 5 per cent of the mass-loss rate. It is assumed that the neutron star accretes matter up to its Eddington limit for helium accretion ($`3\times 10^8`$ $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$); the rest of the transferred matter is lost from the system with the specific orbital angular momentum of the neutron star.

In the limiting case when the mass-transfer rate owing to RLOF is much higher than both the maximum accretion rate and the wind mass-loss rate ($`\dot{M}`$ $`>10^6`$ $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$), we have $`\beta 1`$ and $`\dot{J}_{\mathrm{gwr}}`$ can be neglected. This is generally the case during case BAB and case BB mass transfer discussed below. Then eq. (5) leads to:

$`{\displaystyle \frac{\dot{P}}{P}}={\displaystyle \frac{3q^22q3}{1+q}}{\displaystyle \frac{\dot{M}_{\mathrm{He}}}{M_{\mathrm{He}}}}`$ (11)
From this it can be seen that $`P`$ has a minimum when $`q=1.39`$, i.e. $`M_{\mathrm{He}}1.95`$. For $`q>1.39`$, $`P`$ initially decreases as a result of RLOF, and increases again after $`q<1.39`$. This behaviour can be seen in Sect. 4 and 5.

## 3 The method of calculation

### 3.1 A brief description on the numerical code

We used the Eggleton stellar-evolution code (Eggleton 1971, 1972, 1973). This code uses a self-adaptive, non-Lagrangian mesh-spacing which is a function of local pressure, temperature, Lagrangian mass and radius. It treats both convective and semi-convective mixing as a diffusion process and finds a simultaneous and implicit solution of both the stellar structure equations and the diffusion equations for the chemical composition. The most recent update, including the nuclear reaction network, is described in Pols et al. (1995).

We evolved zero-age main-sequence stars (in Sect. 3.2) assuming a chemical composition of ($`X=0.70,Z=0.02`$) and using a mixing-length parameter of $`\alpha =l/H_p=2.0`$. Convective overshooting is taken into account in the same way as in Pols et al. (1998) using an overshooting constant $`\delta _{\mathrm{ov}}=0.12`$.

The helium stars were evolved assuming a chemical composition of ($`Y=0.98,Z=0.02`$). All calculations were done without enhanced mixing (STD models in Pols 2002). We refer to the latter paper for the evolution of single helium stars related to our calculations.

### 3.2 Initial conditions

#### 3.2.1 Progenitor systems of helium star and neutron star binaries

To obtain estimates of the mass of the helium star and the initial period, we evolved massive main-sequence stars with masses of 8 – 20 $`\mathrm{M}_{\mathrm{}}`$. If the radius of the star is assumed to be equivalent to its Roche radius when RLOF is initiated ($`RR_\mathrm{L}`$), then the separation of the binary system at the onset of the mass-transfer phase, $`a_\mathrm{i}`$, for a given radius $`R`$ can be calculated with eq. (10). With the energy equation for CE evolution (Webbink 1984; de Kool 1990)

$`{\displaystyle \frac{a_\mathrm{f}}{a_\mathrm{i}}}`$ $`=`$ $`{\displaystyle \frac{M_{\mathrm{core}}M_2}{M_{\mathrm{donor}}}}{\displaystyle \frac{1}{M_2+2M_{\mathrm{env}}/(\eta _{\mathrm{CE}}\lambda r_\mathrm{L})}}`$ (12)
we calculated the final separation, $`a_\mathrm{f}`$. Here $`M_{\mathrm{donor}}=M_{\mathrm{core}}+M_{\mathrm{env}}`$, $`M_{\mathrm{core}}`$, and $`M_{\mathrm{env}}`$ are the masses of the donor star, its helium core and envelope, respectively. Due to the stellar wind mass loss (de Jager, Nieuwenhuijzen & van der Hucht 1988; Nieuwenhuijzen & de Jager 1990), $`M_{\mathrm{donor}}`$ might be smaller than $`M_{\mathrm{ZAMS}}`$, the mass of the star on the zero-age main sequence (ZAMS). During the red supergiant phase, the stellar mass is reduced by 3 – 19 per cent (in a 8 $`\mathrm{M}_{\mathrm{}}`$ and a 20 $`\mathrm{M}_{\mathrm{}}`$ stars, respectively). $`M_2`$ is the mass of the companion, which we assumed to be a 1.4 $`\mathrm{M}_{\mathrm{}}`$ neutron star. $`r_\mathrm{L}=R_\mathrm{L}/a_\mathrm{i}`$ is the dimensionless Roche radius of the donor star. We took $`\eta _{\mathrm{CE}}`$, the so-called CE efficiency parameter, to be one. $`\lambda `$, the parameter of the binding energy of the envelope to the core, is approximated as in Dewi & Tauris (2000). We assumed that the binary will survive the CE and spiral-in phase if the final separation is larger than the Roche radius of the helium core, i.e. the helium star does not immediately fill its Roche lobe. Following Taam, King & Ritter’s (2000) argument, we assume that the neutron star does not accrete a significant amount of matter while it is embedded in the donor’s envelope.

The results of the calculations are summarized in Fig. 1, where we have plotted the maximum periods (which are related to the maximum radius) after case B and case C spiral-in. These values are very uncertain, especially for case B, because the parameter $`\lambda `$ used for calculating the post-CE separation depends quite sensitively on the choice of the mass of the helium core. As argued by Tauris & Dewi (2001) approximate lower limits to $`\lambda `$ and to the core mass are obtained by assuming that the core mass is the mass where the abundance of hydrogen is below 0.1. Considering the bottom of convective envelope as the core, results in approximate upper limits of $`\lambda `$ and the core mass. The latter criterion corresponds roughly to the mass where $`X_\mathrm{H}=0.3`$, for the maximum radius during case B from $`M_{\mathrm{ZAMS}}13`$ $`\mathrm{M}_{\mathrm{}}`$. The dotted lines in Fig. 1 mark the maximum periods for case B evolution, calculated with the two extreme $`\lambda `$-values. For case C evolution, only the maximum period (dashed-line) corresponding to the lower limit of $`\lambda `$ is plotted. In stars with $`M_{\mathrm{ZAMS}}14`$ $`\mathrm{M}_{\mathrm{}}`$, helium is ignited before a deep convective envelope has developed, and the radius keeps expanding during core He burning. At the maximum radius for case B (which in this case is not well-defined), the envelope is more strongly bound and the maximum period after case B is correspondingly smaller (see Fig. 1). In this case, the bottom of the convective envelope can no longer be realistically taken as the core mass, so the upper dotted line in Fig. 1 stops at 13 $`\mathrm{M}_{\mathrm{}}`$. Because of these uncertainties, to relax the assumption on the choice of core mass as well as on $`\eta _{\mathrm{CE}}`$ (i.e. the possibility of applying $`\eta _{\mathrm{CE}}`$ other than unity), we did the calculations for a wide range of periods.

#### 3.2.2 Studied types of helium star binary evolution

As Fig. 1 shows, the progenitors of helium star and neutron star binaries were main-sequence star and neutron star systems that evolved according to case B or case C mass transfer. The outcome of CE evolution in case B is an unevolved helium star and neutron star, while that in case C is an evolved helium star (which has developed a carbon core) and neutron star. The systems resulting from case B can be so close that they start RLOF during helium core burning. This is so-called case BA evolution. Wider systems will start RLOF only during helium shell burning around the carbon core. This is so-called case BB evolution. The terminology of the types of evolution used in this paper is summarized in Table 1.

Case C mass transfer will result in evolved helium stars which, after CE evolution, can only start transferring mass after helium core burning is terminated. We are interested in RLOF during helium shell burning, which can be called case CB evolution. Basically, case BB and CB evolutions result in the same type of RLOF from a helium star of about the same mass and period, except in $`M_{\mathrm{He}}4`$ $`\mathrm{M}_{\mathrm{}}`$ where the wind mass loss plays a major role in reducing the mass and changing the structure during helium core burning. We do not separately study case CB evolution in this paper, but assume that it is similar to case BB evolution.

A single helium star expands and reaches the maximum radius during convective carbon shell burning (Habets 1986a; Pols 2002). Although they are not less important, we did not attempt to follow evolutions in which RLOF is initiated during carbon core burning or beyond (case BC or CC evolution). The overall results for case BA mass transfer are shown in Table 2, and those for case BB in Tables 3 and 4.

## 4 Results: Case BA Mass Transfer

Based on the stability of RLOF during helium core burning, we devide case BA mass transfer into three groups, i.e. (i) 1.5 – 2.1 $`\mathrm{M}_{\mathrm{}}`$ helium stars, in which RLOF is stable, (ii) 2.4 – 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium stars, in which RLOF is initiated by a rapid, unstable phase and followed by a slow, stable phase of mass transfer, and (iii) 2.8 – 2.9 $`\mathrm{M}_{\mathrm{}}`$ helium stars, in which the unstable, thermal-timescale RLOF causes the temporary disappearance of the convective core.

In general, case BA mass transfer from helium stars in our calculations resembles case A evolution from main-sequence stars (e.g. Pols 1994; Wellstein, Langer & Braun 2001). In all helium stars of 1.5 – 2.9 $`\mathrm{M}_{\mathrm{}}`$ which undergo case BA mass transfer, the whole envelope is removed from the star, and the remnant, a degenerate CO white dwarf, is never massive enough to ignite carbon.

### 4.1 Roche-lobe overflow from 1.5 – 2.1 $`𝐌_{\mathbf{}}`$ helium stars

During the detached phase the orbit slowly shrinks due to gravitational-wave radiation and the stellar radius slowly increases, until (stable) RLOF starts (point A in Fig. 2). During the first phase of stable, nuclear-timescale RLOF, the star contracts and the energy-production rate in the centre decreases as a result of the loss of mass, which causes the luminosity and effective temperature to decrease. The star remains in thermal equilibrium. The mass of the convective core does not grow as in the case of a single star, but decreases slowly until the core disappears when the helium is exhausted (see Fig. 3). Most of the transferred matter is lost from the system at low specific orbital angular momentum which, although counteracted by gravitational-wave radiation, causes the orbit to slightly expand. Roughly at the same time when the abundance of carbon in the centre is maximum, the period reaches a maximum value, while $`L`$ and $`T_{\mathrm{eff}}`$ reach their minimum values (the reddest point in the HR diagram, point B). After this point, the orbit gradually shrinks. We find that the earlier RLOF is initiated (the smaller the initial period), the more mass is removed from the helium star, the more significant the change in period, and the longer the duration of RLOF.

After the turning point, the star moves to the left in the HR diagram, only an insignificant amount of mass is transferred from the helium star – this is also the case if the orbital period is such that mass transfer starts after the carbon abundance in the centre reaches its maximum. In the latter case, there is no maximum period; the period continuously decreases (while the star shrinks) until the system is detached. The maximum mass-loss rate, the duration of RLOF, and the amount of mass transferred to the neutron star all decrease with decreasing initial mass ratio and with increasing initial period.

As the helium abundance in the core drops below about 0.1, the star contracts inside its Roche lobe, and the system becomes detached (point C). During this detached phase, the orbit is again tightened due to gravitational-wave radiation, but the net change in orbit is insignificant compared to the initial period. The duration of helium core burning is prolonged due to the mass transfer. However, we found that if the initiation of RLOF takes place close to the end of the helium main sequence, helium core burning proceeds marginally faster (by less than 1 per cent) than in a single helium star of the same mass.

After helium is exhausted in the centre and the convective core disappears, helium burning moves to a radiative shell around the CO core. The CO core grows in mass but decreases in radius, which causes the star to expand and move up in the HR diagram. A second phase of stable RLOF, which we call case BAB, then takes place (point D). The central temperature increases. The mass-loss rate, which is higher than that during the helium core burning, increases both with initial mass ratio and with initial period. The orbit expands in all cases due to the effects described in Sect. 2.2. At the end of this mass-transfer phase, when almost the entire envelope is removed from the star, the degeneracy border is crossed. The central temperature decreases as well as the energy-production rate. Although the core is constant in mass, the radius decreases and the star contracts inside its Roche lobe. The system is again detached (point E). The duration of the second RLOF phase decreases with increasing initial period and initial mass ratio.

### 4.2 Roche-lobe overflow from 2.4 – 2.5 $`𝐌_{\mathbf{}}`$ helium stars

RLOF is initially unstable and results in a rapid, thermal-timescale phase of mass transfer. The rapid loss of mass causes a rapid contraction of the star. The mass of the convective core decreases rapidly when the star adjusts to the new total mass. The luminosity and effective temperature decrease, as well as the central temperature and the orbital period, until the mass-loss rate reaches a maximum value. The mass of the convective core then decreases slowly and the star contracts slowly. Mass transfer continues stably on a nuclear timescale. In a 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium star, the transition between the rapid, thermal-timescale phase and the slow, stable mass-transfer phase is marked by a minimum in mass-loss rate (see Fig. 4). The mass of the convective core reaches a minimum value before it increases and then decreases slowly. The further evolution is very similar to that of the lower-mass stars described in Sect. 4.1.

As in the case of 1.5 – 2.1 $`\mathrm{M}_{\mathrm{}}`$ helium stars, the star shrinks inside its Roche lobe before helium is depleted in the centre. The star again fills its Roche lobe when the radius increases after helium is exhausted in the centre (case BAB).

### 4.3 Roche-lobe overflow from 2.8 – 2.9 $`𝐌_{\mathbf{}}`$ helium stars

The mass transfer from 2.8 – 2.9 $`\mathrm{M}_{\mathrm{}}`$ helium stars is also initiated (point 1 in Fig. 2) by a thermal-timescale phase. This time, the rapid mass transfer with $`q2`$ causes the star to lose about 40 per cent of its mass rapidly. The central temperature decreases below the required temperature for helium burning so that the energy production in the centre ceases (point 2, see also Fig. 5). This, in turn, causes rapid contraction of the star and disappearance of the convective helium core. The period decreases rapidly. After the maximum mass-loss rate is reached, the star contracts inside its Roche lobe, the orbit widens due to the loss of mass from the system, and the system is detached (point 3).

After thermal stability is restored (point 4), the convective core reappears, helium is reignited, and the star slowly expands. Wind mass loss is now that of a $``$ 1.7 $`\mathrm{M}_{\mathrm{}}`$ helium star. The orbital period again decreases due to gravitational-wave radiation. In a 2.8 $`\mathrm{M}_{\mathrm{}}`$ helium star in a very close orbit, at the moment when the central carbon abundance reaches a maximum, the stellar radius reaches its Roche radius, and a brief phase of stable RLOF takes place (point 5 to 6). Binaries of the same mass in a slightly wider orbit (as well as those with $`M_{\mathrm{He}}`$ = 2.9 $`\mathrm{M}_{\mathrm{}}`$) do not undergo this second, slow phase of mass transfer since the maximum radius (which relates to the maximum energy-production rate and to the maximum central carbon abundance) is still within its Roche radius. However, the net effect is the same since during the slow phase only an insignificant amount (2 per cent) of mass is removed from the helium star. As in the case of lower-mass stars, the 2.8 – 2.9 $`\mathrm{M}_{\mathrm{}}`$ helium stars also undergo case BAB mass transfer (point 7 to 8) with a lower mass-loss rate than that during the first, rapid RLOF phase.

In more massive helium stars ($`M_{\mathrm{He}}>2.9`$ $`\mathrm{M}_{\mathrm{}}`$), case BA causes the loss of more than half of the stellar mass. The rapid mass loss causes the disappearance of the convective helium core. However, the star is unable to recover its thermal equilibrium. RLOF becomes dynamically unstable and the binary will probably coalesce.

### 4.4 The role of gravitational-wave radiation

Before the star fills its Roche lobe, the orbit is tightened due to gravitational-wave radiation. In order to see the effect of gravitational-wave radiation, we also studied the case in which gravitational-wave emission is turned off. The exclusion of loss of angular momentum due to gravitational-wave radiation will only delay and shorten the first RLOF phase. In this case, during the slow stable RLOF phase, the orbit expands due to the loss of mass from the system, and there is no maximum period at the moment the central carbon abundance reaches its maximum.

The exclusion of gravitational-wave emission in the case of higher-mass helium stars ($`M_{\mathrm{He}}>2.9`$ $`\mathrm{M}_{\mathrm{}}`$) with very close orbits allows the wind mass loss to govern the loss of angular momentum. The orbit is widened such that mass transfer will only take place after the exhaustion of helium in the centre (case BB).

## 5 Results: Case BB Mass Transfer

We devide case BB mass transfer into two mass groups, i.e. (i) 1.5 – 2.1 $`\mathrm{M}_{\mathrm{}}`$ helium stars, in which RLOF is terminated when the degeneracy border is crossed and a heavy, degenerate CO or ONe white dwarf is produced, and (ii) 2.4 – 6.6 $`\mathrm{M}_{\mathrm{}}`$ helium stars, in which an ONe white dwarf or a neutron star is the outcome.

### 5.1 Roche-lobe overflow from 1.5 – 2.1 $`𝐌_{\mathbf{}}`$ helium stars

Case BB mass transfer from 1.5 – 2.1 $`\mathrm{M}_{\mathrm{}}`$ helium stars resembles case BAB in the same range of mass. The exhaustion of helium in the core causes the star to expand and fill its Roche lobe. After the degeneracy boundary is crossed, it contracts rapidly inside its Roche lobe, and the system is detached. The mass-loss rate increases while the duration of RLOF decreases with increasing initial period and initial mass ratio.

After the core becomes degenerate but before the system detaches, carbon burning occurs in the radiative shell with maximum temperature. In 2.1 $`\mathrm{M}_{\mathrm{}}`$ helium stars, off-centre C-ignition takes place in the shell with maximum temperature. At first, carbon burns in a broad radiative shell (and marginally in the centre). When the off-centre convective shell appears, the whole core expands and the region interior to the convective shell cools. C-burning stops completely in the inner core, while convective off-centre burning continues. The calculations break down shortly after this point. It is expected (see Pols 2002 for a discussion) that eventually C-burning reaches the centre and a degenerate ONe core is formed. In the 1.5 and 1.8 $`\mathrm{M}_{\mathrm{}}`$ helium stars, the temperature is never high enough to cause significant carbon burning.

Before RLOF occurs, the gravitational-wave radiation brings the components closer together. In the systems with wider orbits ($`P0\stackrel{d}{.}5`$), the effect of gravitational-wave radiation is weak so that the period remains constant. The evolution of the period is as expected from Sect. 2.2.

We are not always able to follow the evolution until the detachment. However, the drop in mass-loss rate suggests that RLOF is about to end. Almost the entire envelope is removed from the star during RLOF, supporting Delgado & Thomas’ (1981) result. The remnants of 1.5 and 1.8 $`\mathrm{M}_{\mathrm{}}`$ helium stars are degenerate CO white dwarfs with masses of 0.9 – 1.05 $`\mathrm{M}_{\mathrm{}}`$, which are more massive than the remnants of case BAB mass transfer from helium stars in the same range of mass. This trend resembles Pols’ (1994) result for case AB and B mass transfers. Helium stars with initial mass 2.1 $`\mathrm{M}_{\mathrm{}}`$ produce ONe white dwarfs.

### 5.2 Roche-lobe overflow from 2.4 - 6.6 $`𝐌_{\mathbf{}}`$ helium stars

The following classification is based on the time of occurrence of the maximum mass-transfer rate; whether it occurs far before, just before, or during carbon core burning. This can be easily recognized from their tracks in the HR diagram, as are shown in Fig. 6. The close orbit systems are characterized by a small hump (indicated by a star) which marks the transition from the rapid to the slow mass-transfer phase. The wide orbit systems are characterized by the up and down movement (see the right-hand side inset in Fig. 6) which is caused by the decreasing and increasing of the mass-loss rate due to the appearance and disappearance of the C-burning convective shells (except in a 2.4 $`\mathrm{M}_{\mathrm{}}`$ helium star). The intermediate orbit systems show none of the two characteristics.

Before RLOF is initiated, gravitational-wave radiation only plays a role in shrinking the orbit in helium stars with masses of 2.4 – 3.4 $`\mathrm{M}_{\mathrm{}}`$ in close orbits during the early phase of their evolution. In stars of the same masses in wider orbits, gravitational-wave emission is negligible so that the orbit remains constant. Wind mass loss dominates the loss of angular momentum and widens the orbit during the giant phase of the above range of mass, and during the whole evolution of $`M_{\mathrm{He}}3.7`$ $`\mathrm{M}_{\mathrm{}}`$ stars.

#### 5.2.1 Close-orbit systems

In the close-orbit systems, mass transfer takes place initially on the donor star’s thermal timescale. The star contracts and the luminosity and effective temperature – as well as the period – decrease, until the mass-loss rate reaches a maximum, $`\dot{M}_{\mathrm{max},1}`$. $`T_{\mathrm{eff}}`$ then increases until the point which marks the transition between the rapid- and slow-phase of mass transfer (the hump marked by a star in Fig. 6). $`L`$ reaches another minimum, while $`P`$ also decreases to a minimum (point A). The star then moves to the upper left in the HR diagram.

The appearance of a C-burning convective core (or off-centre shell in $`M_{\mathrm{He}}2.5`$ $`\mathrm{M}_{\mathrm{}}`$ stars) causes an increase in mass-loss rate to another maximum, $`\dot{M}_{\mathrm{max},2}`$, which varies from 33 (in a 5.4 $`\mathrm{M}_{\mathrm{}}`$ helium star) to 58 per cent (in a 2.4 $`\mathrm{M}_{\mathrm{}}`$ helium star) of $`\dot{M}_{\mathrm{max},1}`$. After the convective core stops growing in mass in $`M_{\mathrm{He}}4`$ $`\mathrm{M}_{\mathrm{}}`$ stars, during carbon core burning, the mass-loss rate drops, the star moves to the left in the HR diagram and shrinks inside its Roche lobe (point B, see the left-hand side inset in Fig. 6). The period is relatively constant during this detached phase.

When the carbon abundance in the convective core drops below 0.1, the star expands again, moves to the right in HR diagram, and the orbit expands.222We could not follow the subsequent evolution in $`M_{\mathrm{He}}2.5`$ $`\mathrm{M}_{\mathrm{}}`$ stars. However, we expect that these stars will undergo the same phenomenon as is found in the more massive ones. The second phase of RLOF occurs after the C-burning convective core has vanished (point C). The mass-loss rate increases rapidly and reaches a maximum, $`\dot{M}_{\mathrm{max},3}`$ (point D), when the first C-burning convective shell appears. When this shell disappears, the mass-loss rate drops again.

In $`M_{\mathrm{He}}2.9`$ $`\mathrm{M}_{\mathrm{}}`$ stars, the radius decreases and the components are again detached (point E). The star moves to the left again at constant period until it reaches a minimum radius when the second C-burning convective shell appears. The star fills its Roche lobe for the third time (point F) after the second convective shell has vanished. $`\dot{M}`$ increases and reaches a maximum value, $`\dot{M}_{\mathrm{max},4}`$. The second RLOF lasts shorter than the first one, and the third one shorter than the second, because the evolution of the interior region proceeds faster (Habets 1986a).

In 3.1 – 4 $`\mathrm{M}_{\mathrm{}}`$ helium stars, RLOF takes place in two episodes. After the first C-burning convective shell disappears, the mass-loss rate decreases to a minimum, but RLOF continues without a detached phase. In 3.6 – 4 $`\mathrm{M}_{\mathrm{}}`$ helium stars, the first shrinkage of the star inside its Roche lobe takes place later than in the lower-mass stars, i.e. when carbon is being depleted in the core. The radius reaches its minimum value after the core has vanished. With the appearance of the first C-burning convective shell, the star fills its Roche lobe for the second time.

In $`M_{\mathrm{He}}>4`$ $`\mathrm{M}_{\mathrm{}}`$ stars, $`\dot{M}`$ also decreases after reaching $`\dot{M}_{\mathrm{max},2}`$, until carbon is exhausted in the core. The behaviour of $`\dot{M}`$ subsequent to carbon core burning is erratic, but it never exceeds $`\dot{M}_{\mathrm{max},1}`$. In this range of mass, the star never shrinks inside its Roche lobe.

Why detached phases occur in one case and not in the other, and why the detachments take place in different phases of carbon burning for different cases, can be understood as follows. The start of each convective carbon burning phase (both in the core and in subsequent shells) causes the core to expand. As a result, the envelope tends to contract and the star tends to shrink inside its Roche lobe. In between the convective phases, the core contracts and the envelope expands again. However, the reaction of the stellar surface is delayed by the thermal diffusion timescale between the core and the surface. This timescale is roughly constant for each star, and does not vary strongly with mass. On the other hand, the timescale for C-burning decreases strongly with each successive convective phase, as well as with increasing mass, and eventually become much shorter than the thermal timescale. As a result, the detachments are delayed or even absent for advanced phases and/or large masses.

The transition from rapid to slow mass transfer is not found in $`M_{\mathrm{He}}>5.4`$ $`\mathrm{M}_{\mathrm{}}`$ stars. Stars in this mass range with periods shorter than those in Table 4 go through dynamically unstable mass transfer.

#### 5.2.2 Intermediate-orbit systems

In the intermediate-orbit systems, RLOF occurs on the thermal timescale, the orbit shrinks, and the maximum mass-transfer rate, $`\dot{M}_{\mathrm{max},1}`$, is reached not long before the appearance of the C-burning convective core (or off-centre convective shell in 2.4 – 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium stars). The star moves afterwards to the upper left in the HR diagram (see Fig. 6). Unlike in the close-orbit systems, the mass-transfer rate does not increase when the convective core appears. Along with the increase in $`L`$ and $`T_{\mathrm{eff}}`$, $`\dot{M}`$ decreases.

In 2.4 – 2.9 $`\mathrm{M}_{\mathrm{}}`$ helium stars, $`\dot{M}`$ decreases until the star shrinks inside its Roche lobe during carbon core (or off-centre shell) burning, and the system is detached (point a, see the middle inset of Fig. 6 and Fig. 8). The period increases shortly after $`\dot{M}_{\mathrm{max},1}`$ is reached, and is relatively constant during the detached phase. When the carbon abundance in the core (or the off-centre shell) drops below 0.1, the star reaches the bluest point (In 2.4 – 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium stars, at this point the off-centre shell penetrates to the centre). The star expands, moves back to the right in the HR diagram, and fills its Roche lobe after the C-burning convective core has disappeared (point b). During the second RLOF, $`T_{\mathrm{eff}}`$ and $`L`$ decrease while the period increases. After the first C-burning convective shell disappears, the mass-transfer rate reaches a maximum, $`\dot{M}_{\mathrm{max},2}`$ (point c), before decreasing again (see Fig. 8).

In 3.1 – 5.8 $`\mathrm{M}_{\mathrm{}}`$ helium stars in intermediate orbits, RLOF takes place without detached stages. After $`\dot{M}_{\mathrm{max},1}`$, $`\dot{M}`$ goes through a minimum and then increases again. Also due to the difference in thermal timescale between the core and the envelope, the occurence of the minimum $`\dot{M}`$ (particularly in this range of mass) varies. Like in the close-orbit systems, the period decreases throughout the whole mass-transfer phase.

#### 5.2.3 Wide-orbit systems

Thermal-timescale RLOF is most rapid in wide-orbit systems, and the maximum mass-loss rate, $`\dot{M}_{\mathrm{max},1}`$, is reached during the carbon core (or off-centre shell) burning. The mass-transfer rate then decreases and reaches a minimum value (in a 2.4 $`\mathrm{M}_{\mathrm{}}`$ helium star, the decrease in $`\dot{M}`$ leads to a detachment), before reaching another maximum, $`\dot{M}_{\mathrm{max},2}`$. The increases and decreases in $`\dot{M}`$ can be seen in the HR diagram as the up and down movements (see Fig. 6). Except in a 2.4 $`\mathrm{M}_{\mathrm{}}`$ helium star, RLOF in wide-orbit systems takes place without any detached stages. Again, only in $`M_{\mathrm{He}}<2.9`$ $`\mathrm{M}_{\mathrm{}}`$ stars does the period reach a minimum value shortly after $`\dot{M}_{\mathrm{max},1}`$.

#### 5.2.4 The He-burning convective envelope and common-envelope evolution

As can be seen in the right-hand panels of Fig. 9, a He-burning convective shell appears in the envelope before the appearance of the third C-burning convective shell. Such a shell also appears in previous helium star calculations (e.g. Habets 1986a, 1986b; Nomoto 1984, 1987; Pols 2002). Because of the small envelope mass, this shell almost penetrates to the surface after the disappearance of the third C-burning convective shell. In response to mass loss from a convective envelope, the star expands, $`\dot{M}`$ increases enormously by 2 orders of magnitude, and the orbit widens. $`L`$ and $`T_{\mathrm{eff}}`$ decrease and the star moves to the lower right in the HR diagram (see Fig. 6). In a 3.8 $`\mathrm{M}_{\mathrm{}}`$ helium star, the He-burning convective shell appears and penetrates to the surface before the second C-burning convective shell appears.

We are not always able to follow the evolution of the helium stars beyond the appearance of the second C-burning convective shell. However, we found that the behaviour occurs for those masses and orbits that we were able to compute the furthest (see the c notation in Table 4) in helium stars less massive than 4 $`\mathrm{M}_{\mathrm{}}`$. Therefore, we expect that this is a typical phenomenon of the late stage of a case BB mass transfer from $`M_{\mathrm{He}}<4`$ $`\mathrm{M}_{\mathrm{}}`$ stars (in RLOF initiated during carbon shell burning (Habets 1986a), the envelope is thick so that the expansion of He-burning convective shell occurs deep inside the envelope). In $`M_{\mathrm{He}}4`$ $`\mathrm{M}_{\mathrm{}}`$ the start of a convective shell is only found in close orbits, but not in the wide orbits athough the calculations were done up to the same stage of evolution. However, we cannot exclude that such a convective shell does not eventually develop in these more massive stars as well.

The outcome of mass transfer in the presence of a convective helium envelope is not clear. Mass transfer might become dynamically unstable which would lead to a CE and spiral-in phase. Although we do not find a convective envelope in helium stars more massive than 4 $`\mathrm{M}_{\mathrm{}}`$ in wider orbits, the large mass ratio ($`q2.6`$ at the onset of RLOF) and the high mass-loss rate ($`\dot{M}`$ $`>10^4`$ $`\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$) suggest that a CE and spiral-in phase will very likely also take place. However, the high luminosity due to neutrino emission ($`L_\nu >10^7`$ $`\mathrm{L}_{\mathrm{}}`$) suggests that the core will collapse within a few decades after the end of the calculation. The closeness to collapse is supported by the fact that Ne-burning has already started radiatively off-centre when the fourth C-burning convective shell appears. The explosion will probably be observed as a type Ic supernova. Nevertheless, should the two components undergo a spiral-in, a merger will not necessarily occur. A crude calculation of the gravitational binding energy of the convective envelope to the CO-core boundary of a 3.2 $`\mathrm{M}_{\mathrm{}}`$ helium star when the penetration takes place yields $`\lambda `$ = 0.123. With $`\eta _{\mathrm{CE}}=1`$, $`M_{\mathrm{He}}`$ = 1.98 $`\mathrm{M}_{\mathrm{}}`$, $`M_{\mathrm{core}}`$ = 1.62 $`\mathrm{M}_{\mathrm{}}`$ and $`P_{\mathrm{orb}}=0\stackrel{d}{.}22`$, we obtain from eq. (12) a final separation of 0.0729 $`\mathrm{R}_{\mathrm{}}`$ (period of $`1\stackrel{m}{.}89`$), which is still larger than the radius of the core (0.0275 $`\mathrm{R}_{\mathrm{}}`$). Without helium left in the envelope, the explosion will be observed as a type Ic supernova. The close post-SN separation subjects the two neutron stars to a strong loss of gravitational-wave radiation, causing the two compact stars to merge on a timescale of 2 500 years, which may produce a gamma-ray burst.

#### 5.2.5 Final masses and outcome

The final masses and periods are shown in Table 4. Here final corresponds to the last stage of our calculations, which varies from the appearance of the first C-burning convective shell to the penetration of the He-burning convective shell to the surface. We found that after the carbon convective core has disappeared, only a small amount of mass is transferred to the neutron star. The period alters insignificantly. Hence, the values in Table 4 can be considered as the final masses and periods prior to SN explosion, unless a spiral-in phase occurs as discussed in Sect. 5.2.4.

With increasing initial period and decreasing initial mass, the amount of transferred matter decreases. In 2.4 – 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium stars, the mass of the remnant ranges from 1.3 – 1.6 $`\mathrm{M}_{\mathrm{}}`$, but the mass of the CO core (1.21 – 1.34 $`\mathrm{M}_{\mathrm{}}`$) never exceeds the Chandrasekhar limiting mass, $`M_{\mathrm{Ch}}`$. These stars will develop degenerate ONe cores. The future of these systems is not clear since they could not be computed through carbon burning. If the ONe core can grow to 1.37 – 1.39 $`\mathrm{M}_{\mathrm{}}`$, they will undergo a SN explosion triggered by electron capture on Ne and Mg (e.g. Nomoto 1984, 1987), otherwise they will become heavy ONe white dwarfs.

Stars with $`M_{\mathrm{He}}2.62.8`$ $`\mathrm{M}_{\mathrm{}}`$ which undergo case BB mass transfer develop CO cores more massive than $`M_{\mathrm{Ch}}`$ (where the actual lower mass limit probably depends on the orbital period). These stars will eventually complete neon, oxygen and silicon burning and form iron cores. The collapse of the Fe core by photodisintegration triggers a SN explosion, leaving a neutron-star remnant.

The final amount of helium left in the envelope probably determines whether the explosion will be a Type Ib or a Type Ic supernova. The main observational criterion to distinguish between these types is the detection of helium in the spectrum. Helium is generally absent in SNe Ic, although a small amount may be present (Filippenko et al. 1995). The maximum allowed amount of He is very uncertain, however, and may also depend on other factors such as mixing during the explosion (Woosley & Eastman 1997).

For low-mass He stars and/or close orbits, we find there may be as little as 0.04 $`\mathrm{M}_{\mathrm{}}`$ of He left, and generally $`<`$ 0.10 $`\mathrm{M}_{\mathrm{}}`$ (see Table 4). This is less than is achieved by stellar-wind mass loss in massive He-star models (Woosley, Langer & Weaver 1995; Wellstein & Langer 1999; Pols 2002) and also less than obtained from conservative case ABB or case BB mass transfer models by Wellstein & Langer (1999). We conclude that case BB mass transfer to a compact companion in a very-close binary is the most efficient way to produce almost bare CO cores prior to explosion. These stars will almost certainly produce a Type Ic supernova. Such a model was first suggested by Nomoto et al. (1994) as the progenitor of the Type Ic SN 1994I, and shown to match the observed lightcurve (Iwamoto et al. 1994). The amount of He left in the envelope increases with increasing initial mass and initial period. Probably, the more massive case BB remnants and those from the widest orbits undergo a Type Ib explosion, as do helium stars that avoid case BB mass transfer (Pols 2002).

### 5.3 The binary millisecond pulsars and double neutron-star pulsars

The outcome of case BB mass transfer from 1.5 – 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium stars are heavy white dwarfs with a neutron-star companion. Their positions in the period-mass plane (solid- and open-circles in Fig. 10) cover the region where the existence of binary millisecond pulsars with heavy white-dwarf companions cannot be explained by extreme mass transfer (Tauris, van den Heuvel & Savonije 2000). Systems like J1756-5322 ($`P_{\mathrm{orb}}=0\stackrel{d}{.}453`$, $`M_{\mathrm{WD}}=0.683`$ $`\mathrm{M}_{\mathrm{}}`$ with the assumption that $`i=60^{}`$) are not produced by our nor Tauris et al.’s calculations. Therefore, they must have come from stars with $`5<M_{\mathrm{ZAMS}}/\mathrm{M}_{\mathrm{}}<8`$.

Case BB mass transfer from $`M_{\mathrm{He}}2.62.8`$ $`\mathrm{M}_{\mathrm{}}`$ produces neutron-star binaries. Attempts to find the progenitors of the three known galactic double neutron-star binaries (outside the globular cluster) have been carried out e.g. by Yamaoka, Shigeyama & Nomoto (1993), Fryer & Kalogera (1997), Bagot (1997), and recently Francischelli, Wijers & Brown (2001). Their calculations were approximated by the evolution of single helium stars, and it is often assumed that the helium star does not fill its Roche lobe (i.e. $`R_{\mathrm{He}}<R_{\mathrm{L},\mathrm{He}}`$) in order to avoid an unstable phase of mass transfer. The orbital parameters of the observed binary pulsars are presented in Table 5. Here we will try to see if our helium star calculations can explain their existence.

If we simply assume that the radius of the pre-SN orbit (the helium-star binary), $`a_\mathrm{o}`$, must be inside the radius turning point of the post-SN orbit, $`a_\mathrm{p}`$, (Flannery & van den Heuvel 1975) i.e.,

$`(1e)a_\mathrm{p}a_\mathrm{o}(1+e)a_\mathrm{p}`$ (13)
where $`e`$ is the eccentricity of the binary pulsar, then from Kepler’s third law we have the following condition

$`(1e)^3M_\mathrm{p}P_\mathrm{p}^2M_\mathrm{o}P_\mathrm{o}^2(1+e)^3M_\mathrm{p}P_\mathrm{p}^2`$ (14)
where $`M_\mathrm{o}`$, $`M_\mathrm{p}`$ are the total pre- and post-SN masses; $`P_\mathrm{o}`$, $`P_\mathrm{p}`$ are the pre- and post-SN periods.

We plotted the allowed minimum and maximum pre-SN periods according to eq. (14) for each double neutron-star binary in Fig. 10, under the assumption that the pulsar companion is a 1.4 $`\mathrm{M}_{\mathrm{}}`$ neutron star. Fig. 10 shows that systems like B1913+16 are easy to produce from systems with initial masses of 2.4 – 6.7 $`\mathrm{M}_{\mathrm{}}`$ in initial periods of $`0\stackrel{d}{.}1`$$`0\stackrel{d}{.}7`$, while systems like B1534+12 can be made from helium stars with the same range of mass in a narrower range of period ($`0\stackrel{d}{.}3`$$`0\stackrel{d}{.}6`$). None of the case BB models we have computed can produce B1518+49-like systems. The remnant of the 2.4 $`\mathrm{M}_{\mathrm{}}`$ system (indicated by the arrow in Fig. 10) is a ONe white dwarf/neutron star binary in a circular orbit. A case BC evolution might be able to explain the formation of this system, or a $`M_{\mathrm{He}}3`$ $`\mathrm{M}_{\mathrm{}}`$ star, which avoids RLOF altogether (the shaded area in Fig. 10).

The region to the left of the solid line in Fig. 10 produces white dwarfs. The area below the solid curve will undergo case BA or dynamically unstable mass transfer (see Sect. 4). Therefore, these regions can be considered as the zone of avoidance for the formation of double neutron stars.

## 6 Conclusions

We have done calculations of helium stars with masses in the interval of 1.5 – 6.7 $`\mathrm{M}_{\mathrm{}}`$ with a 1.4 $`\mathrm{M}_{\mathrm{}}`$ neutron star companion. We allowed the helium star to fill its Roche lobe during helium core burning (case BA) and helium shell burning (case BB mass transfer). The final fate of the helium star as a function of its initial mass and period is summarized in Fig. 11.

Case BA mass transfer from helium stars with masses of 1.5 – 2.9 $`\mathrm{M}_{\mathrm{}}`$ occurs during helium core burning and stops when the central abundance of helium drops below 0.1 and the star contracts. In helium stars less massive than 2.4 $`\mathrm{M}_{\mathrm{}}`$, RLOF is stable and takes place on the nuclear timescale. For $`2.4M_{\mathrm{He}}/\mathrm{M}_{\mathrm{}}2.9`$, a rapid, thermal-timescale phase of RLOF is followed by a stable, slow phase of mass transfer on the nuclear timescale. Driven by radius expansion, a stable phase of case BAB mass transfer occurs during helium shell burning. This phase of mass transfer stops when the degeneracy border is crossed and the star contracts. During RLOF, the entire envelope is removed from the helium star. The remnants of case BA mass transfer are heavy, degenerate CO white dwarfs. This implies that a star with initial mass as large as 12 $`\mathrm{M}_{\mathrm{}}`$ still can become a white dwarf, as also found by e.g. Wellstein et al. (2001) in stars that go through conservative contact-free case A evolution.

Stable case BB mass transfer from 1.5 – 2.1 $`\mathrm{M}_{\mathrm{}}`$ stars takes place on the nuclear timescale and stops when the degeneracy border is crossed and the star contracts. The entire envelope is removed during RLOF, and heavy, CO (from 1.5 – 1.8 $`\mathrm{M}_{\mathrm{}}`$) or ONe (from 2.1 $`\mathrm{M}_{\mathrm{}}`$ helium stars) white dwarfs are the remnants in this range of mass. Case BB mass transfer from 2.4 – 2.5 $`\mathrm{M}_{\mathrm{}}`$ helium stars produces ONe white dwarfs and that from more massive ones produces neutron stars. RLOF from $`M_{\mathrm{He}}2.4`$ $`\mathrm{M}_{\mathrm{}}`$ takes place on the thermal timescale. In close-orbit systems, this is followed by a slow, stable phase of mass transfer. There is a tendency that the helium envelope becomes convective in helium stars less massive than 4 $`\mathrm{M}_{\mathrm{}}`$. Mass transfer might become dynamically unstable in this case. However, the star probably collapses before a spiral-in can occur. Even if a CE and spiral-in phase does occur, the system is likely to survive. In both cases, a double neutron-star binary is formed.

From the amount of helium left in the envelope, we suggest that helium stars of high mass or in a wide orbit produce type Ib SNe. Type Ic SNe come from helium stars of lower mass or in a close orbit. We are also able to provide a zone of avoidance for the formation of double neutron stars.

## Acknowledgements

This work was sponsored by NWO Spinoza Grant 08-0 to E. P. J. van den Heuvel. It is a pleasure to thank Norbert Langer, Gijs Nelemans, Winardi Sutantyo, Thomas Tauris, and Lev Yungelson whose valuable comments and suggestions improve the early version of the manuscript; and Philipp Podsiadlowski for his constructive referee comments.


format: text

text:
# The multiperiodicity of the 𝛾 Doradus stars HD 224945 and HD 224638 as detected from a multisite campaignBased on observations partially collected at ESO-La Silla (Proposals 54.E-018 and 56.E-0308)

## 1 Introduction

The variability of HD 224638$``$BT Psc ($`V`$=7.5, F1 V) and HD 224945$``$BU Psc ($`V`$=6.93, F0 V) was announced by Mantegazza & Poretti (1991), as a by–product of the monitoring of the $`\delta `$ Sct star HD 224639$``$BH Psc. Both stars had been used as comparison stars in the first observing run devoted to BH Psc. They increased the number of known $`F`$–type stars located close to the low–temperature edge of the Cepheid instability strip which exhibit small amplitude variability on time scales of several hours, usually longer than the length of a night of observation and therefore easily detectable only when used as comparison stars for short–period variables. At that time, the debate on the nature of these variations was divided between spot activity (observed periods as rotational periods) and pulsation (nonradial $`g`$–modes). Mantegazza et al. (1994, hereinafter Paper I) tried to explain the complicated light behaviour of HD 224639 and HD 224945 in the simplest way possible, by means of periodicities shaped as double– or triple–wave curves. Balona et al. (1994) reported on the multiperiodicity of $`\gamma `$ Dor, giving decisive evidence in favour of pulsation. Also, as a result of combined photometry and line profile variations for 9 Aurigae (Krisciunas et al. 1995) and $`\gamma `$ Dor (Balona et al. 1996), plus the preliminary results on HD 224945 (Poretti et al. 1996), the hypothesis of variability caused by $`g`$–modes achieved a wide consensus. The properties of this new class of variable stars were delineated step–by–step through observational efforts and have been summarized by Kaye et al. (1999), while the problem of the driving mechanism and the excitation of $`g`$–modes constitutes the target of continuing theoretical investigations. As a result, $`\gamma `$ Dor variables are now considered intermediate between $`AF`$ pulsators and solar–type stars, finding a special place in the main programmes of asteroseismologic missions such as corot and mons.

On the basis of this progress, it is worth investigating in more detail the pulsational behaviour of HD 224639 and HD 224945. Here we present the results of a multisite campaign carried out in October 1995.

## 2 Observations

HD 224945 and HD 224638 are located very close to the celestial equator and therefore can be monitored from both the northern and southern hemispheres. As comparison stars we used the same two stars as in the 1991 campaign (HD 225086 and HD 200), since they proved to be stable within a few mmag. Five telescopes were used for this multisite campaign:

1. The European Southern Observatory 50–cm telescope located in La Silla (Chile), equipped with a photon–counting photometer (EMI 9789 QB photomultiplier) and $`B`$ and $`V`$ filters. The observer was E. Poretti.
2. The 48–cm telescope, located at Ege University Observatory, was equipped with a solid–state photometer and $`B`$ and $`V`$ filters. The observer team was led by M. C. Akan.
3. The 50–cm telescope located in Sutherland, South African Astronomical Observatory, was equipped with a photon–counting photometer and $`B`$ and $`V`$ filters. The observer was C. Koen.
4. The 61–cm telescope located on the Mauna Kea (Hawaii) was equipped with a photon–counting photometer and standard $`B`$ and $`V`$ filters. The observers were K. Krisciunas, R. Crowe and M. Wilcox.
5. The 90–cm telescope located at Sierra Nevada Observatory (Spain). The observers were E. Rodríguez and S. Martín. Photometry was performed in the $`uvby`$ system, but only $`v`$ and $`y`$ measurements are discussed here, as the more compatible with the $`B`$ and $`V`$ ones, respectively.

Earlier observations were performed in 1994 at European Southern Observatory, using the ESO 50–cm telescope equipped with the same instrumentation used in 1995. On that occasion the observer was M. Bossi; the results obtained on BH Psc are reported by Mantegazza et al. (1996). The measurements were distributed over 13 nights, with a total useful observing time of 64 hours and a baseline of 14.9 d.

## 3 Data reduction and analysis

The original measurements were provisionally reduced by the observers themselves to check the data quality and for comparison purposes. The original raw data were then re–reduced by means of the same algorithm; this homogenous procedure allowed us to avoid the introduction of local drifts caused by routines using different interpolation formulae and/or different methods of applying extinction corrections. Since the weather was not particularly favourable during the campaign, the rejection of some parts of nights characterized by unfavourable weather conditions was decided on the basis of a uniform criterion. We also decided to use the code which allows us to calculate instantaneous values of the extinction coefficients (Poretti & Zerbi 1993) since in some cases changes were expected to occur.

Table 1 lists some details regarding observing runs at each site. As can be noted, the mean values of the magnitude differences between the two comparison stars show the presence of systematic shifts (up to 0.03 mag). This is not surprising, taking into account that different photomultipliers and filters were used. The systematic shifts between HD 225086 and HD 200 cannot be used to correct magnitude differences between HD 225086 and the two program stars since they have different colours from HD 200; uncertainties of the colour transformations are larger than the precision we need. Hence, to apply this correction we exploited the very useful circumstance that the spread in longitude between different observing sites is not larger than the length of the night. There were always some overlapping segments of the light curve, which allowed us to determine the amounts of the corrections in a straightforward way.

It should also be noted that the Turkish measurements show a scatter larger than all the others. At first, all the Turkish data were considered, but in the last step of the analysis (once it was established that rapid variations are not present in the light curves) those measurements were averaged in groups of 4–5. Thus, the final set of data contains the original measurements from all the observatories except for the averaged values from Turkey. The standard deviations of the $`\mathrm{\Delta }B`$ and $`\mathrm{\Delta }V`$ values between the comparison stars are 3.7 and 3.8 mmag, respectively.

Figure 1 shows a part of the $`B`$ light curves of HD 224945 and HD 224638. The fits derived in the next sections are also shown. The original measurements can be requested from the authors.

Since data in the final form circulated amongst all the participants, several co–authors analysed the time series independently and by different period search algorithms (least–squares, clean, dft, …) and it was particularly satisfying to see that they detected the same terms, even if, owing to the complexity of the light variations, the terms did not always stand out in a clear way.

We present here the analysis carried out by using the least–squares method (Vaniĉek 1971) used by the Merate group in the analysis of $`\delta `$ Sct star light curves. This method has the advantage of not using any data prewhitening since only the frequency values previously found are considered as input values (known constituents; k.c.’s); their amplitude and phase values are recalculated as unknowns in the new searches. That means that after the detection of the $`f_1`$ term, only the frequency value $`f_1`$was considered as established (i.e. a k.c.) and in the second search the unknowns were $`\mathrm{\Delta }m_o,A_1,\varphi _1,f_2,A_2,\varphi _2`$. The ordinates of the power spectra show the reduction factor

$$\mathrm{Red}.\mathrm{Factor}=1\frac{\sigma _{\mathrm{fin}}^2}{\sigma _{\mathrm{in}}^2}$$
(1)
Moreover, we can present the term detection step–by–step. The frequency values were refined after each new detection. We can also obtain useful evidence concerning the amplitude by applying a least–squares fit to the datasets. The interpolating formula we used is

$$\mathrm{\Delta }m(t)=\mathrm{\Delta }m_o+\underset{i}{\overset{N}{}}A_i\mathrm{cos}[2\pi f_i(tT_o)+\varphi _i]$$
(2)
Since the 1995 observations are concentrated in an interval of 15 d only, the frequency resolution is 1.5/$`\mathrm{\Delta }`$T=0.10 cd-1. However, the absence of relevant aliases makes the detection of the excited terms very simple. We note that the frequencies mentioned in the text and figures are given to better accuracies that the formal frequency resolution, because they were determined by a least–squares procedure. The formal least–squares standard errors on the frequencies no doubt understimate the true uncertainties – aside from problems pointed out by e.g. Montgomery & O’Donoghue (1999), it is clear that we have not always resolved close frequencies, and that some unidentified signals may remain in the data.

## 4 Frequency analysis of the HD 224945 data

### 4.1 The 1995 campaign

We collected 442 $`V`$ and 417 $`B`$ measurements of HD 224945, considering the Turkish measurements as binned points. The results of the frequency analysis of the $`B`$ data are shown in Fig. 2. In each panel the horizontal line indicates the level for $`S/N=4.0`$, i.e., the limit usually accepted for significance (Kuschnig et al. 1997).

The top panel shows the spectrum obtained without any k.c.: the peak at 3.00 cd-1 stands out clearly, and the alias structure mimics very well the spectral window (see also Fig. 1 in Poretti et al. 1996). In the discussion of the 1991 data (Paper I), the alias at 2.00 cd-1 was erroneously preferred. The reality of this term close to a multiple of 1 cd-1 was discussed in Paper I. We note that it is also detected in the new campaign data for HD 224945, but in neither the HD 224638, nor the HD 200, time series; its physical presence in the HD 224945 data is therefore certain. When introducing the 3.00 cd-1 frequency as k.c., the power spectrum shows an almost flat pattern, indicating that the amplitude of the remaining terms is smaller. However, two peaks reach the acceptance level: the first is located at 1.16 cd-1 (second panel), the second at 2.84 cd-1 (third panel). After introducing these three terms as k.c.’s, the detection of further terms becomes delicate. The highest peak in the fourth panel is at 2.42 cd-1 and after that the power spectrum does not show a really dominant peak (bottom panel).

The analysis of the $`V`$ measurements yields the same results, except with a slight difference in the frequency of the 2.84 cd-1 term. A value close to 2.77 cd-1 is preferred. Since the peak at 2.84 cd-1 is broad (see third panel of Fig. 2), the presence of another undetectable term combined with the noise distribution could be responsible for the difference between the two values. In the least–squares solution we will consider the 2.84 cd-1 term as the only well–established term.

These independent terms describe a different scenario from that in Paper I. In that case the single–site observations make the power spectrum very complicated, with a great uncertainty between a peak and its $`\pm `$1 cd-1 aliases. This uncertainty led us to propose a solution based on two periodicities each having a triple–wave shape, as it was considered the simplest. Thanks to the multisite observations, we now know that many pulsational modes are simultaneously excited.

### 4.2 The 1994 and 1991 observations

The frequency analysis of the $`B`$ data collected in the 1994 season is much more complicated than that of the multisite campaign. In spite of that, we could detect the 3.00 and the 2.84 cd-1 terms. When introducing these two terms as k.c.’s, the highest peak is at 2.27 cd-1. This peak is significantly different than the 2.42 cd-1 detected in the 1995 multisite campaign. The search for a fourth component is not easy since we get numerous peaks: however, we could recognize the 1.16 cd-1 term and its aliases. We considered it as a real mode on the basis of the results of the 1995 campaign. After that, the new search showed evidence of a complicated structure centered at 1.66 cd-1. We have no way of selecting in a reliable way the true term amongst this peak and the aliased ones. Moreover, the S/N of this term is less than 3.0. The $`V`$ data allowed us to identify the same terms, supporting the presence of the new term at 2.27 cd-1 .

In the 1991 dataset the detection order of the terms is 3.00, 2.84 and 2.31 cd-1 : the only difference with respect to the previous least–squares search (see the beginning of § 3.2 in Paper I) is the identification of the first term as 3.00 cd-1 instead of 2.00 cd-1 . No other term can be detected after considering this triplet as k.c.’s; in particular there is no evidence of the presence of the 1.16 cd-1 term.

It should be emphasized that it would be very difficult to solve the 1994 and 1991 light curves without having the way shown by the results of the 1995 campaign. In turn, the similar structures observed in the three independent datasets corroborate the proposed identification of the excited frequencies, even if some of them are at the limit of the acceptance level.

### 4.3 Frequency refinement and least–squares fitting

The much longer time interval covered by the 1991 data (the two runs span 14 and 10 days, respectively, and they are separated by a gap of 10 days) allowed us to use these data to refine the frequency values. First, we calculated the best solution using all the 1991 data. Then we calculated the amplitudes for each 1991 subset keeping fixed frequencies and phases.

Before applying these frequency values to other datasets, we verified the correspondence between some terms. We considered the 2.27 cd-1 term detected in the 1994 data as the same detected at 2.31 cd-1 in the 1991 campaign. However, we note that the least–squares solution yields a better fit when considering different values for different seasons. On the other hand, we cannot consider the 2.42 cd-1 term detected in 1995 as coincident with the 2.31 cd-1 term. We suggest that there are many terms excited (simultaneously or not) in the narrow range 2.2–2.5 cd-1and their detection is not an easy task. We considered the following five terms as independent modes detected in the three campaigns: 3.003, 2.836, 2.313, 2.424 and 1.160 cd-1. We also note the almost 2:1 ratio between the 2.313 and 1.160 cd-1 terms: the similar amplitude does not support a monoperiodic contribution to the light curve as it should be very asymmetric. The resonance effect is more plausible, even if it should be noted that the 1.160 cd-1 term was not observed in the 1991 dataset.

Using the four frequencies above, we calculated the amplitudes from the 1994 and 1995 datasets. Tab. 2 summarizes the results. The formal errors of the amplitudes are about 0.4 mmag. Looking at the differences between the amplitudes of the two 1991 subsets we can see that they are marginally significant, being in the interval 0.7–0.9 mmag. This very close similarity demonstrates that solutions obtained from short runs are self–consistent. The dense time coverage in each dataset ensures the reliability of the amplitude determination. Therefore, we can look at the differences in the amplitudes between the 1991, 1994 and 1995 datasets with greater confidence.

We also analyzed in frequency the data obtained by merging the 1991, 1994 and 1995 datasets. For each detected term, it is not possible to select the true value owing to the numerous aliases separated by integer values of $`\pm `$1 cy-1; therefore, in the case of HD 224945 we cannot propose a full, comprehensive solution and evaluate in a reliable way the phase coherence of the pulsation. The amplitude variability is a further complication. It should be noted that an attempt at such an analysis detected two peaks at 2.27 and 2.31 cd-1, suggesting that these two peaks are related to two independent terms (see the above comparison of the 1994 and 1991 results).

The sum of the squared amplitudes are the same in 1991 and 1994 (61 and 58 mmag2, respectively), while it is lower in 1995 (50 mmag2). However, the amplitude of a mode can change dramatically: the 2.836 cd-1 term halved its $`B`$ amplitude from 1991 to 1995 while the 2.313 cd-1 term disappeared. As a matter of fact, in the 1991 dataset there are three terms having the same amplitude, while in the 1995 ones the 3.003 cd-1 term is dominating. The amplitude variability looks like a well–established fact for this $`\gamma `$ Dor star.

## 5 Frequency analysis of the HD 224638 data

### 5.1 The 1995 campaign

The analysis of the most recent $`B`$ and $`V`$ measurements of HD 224638 identified the same terms as before, i.e., 1.627, 1.368, 1.697, 1.565 and 1.145 cd-1 (Fig. 3). Also in this case the multisite observations allowed us to identify the terms without any ambiguity. Differently than the case of HD 224945, there are several very close frequencies: 1.627 (top panel), 1.697 (third panel), and 1.565 cd-1 (fourth panel). The separations between the two side peaks and the central one are 0.07 and 0.06 cd-1; these values are comparable to the half width at the half maximum, i.e. 1/$`\mathrm{\Delta }T`$=0.07 cd-1. Indeed, in the power spectra the peaks are not well resolved, but they appear to be double and/or enlarged. In particular, the double peaks appearing in the third and fourth panels of Fig. 3 should be noted. This reveals the presence of overlapping peaks. Two peaks flanking a central one at the limit of the frequency resolution can be generated by a single term modulated in amplitude; however, the amplitudes of the side peaks are comparable to that of the central one and this fact supports the reality of the triplet. The analysis of the 1991 dataset definitely confirms the hypothesis of the three independent modes (see next subsection). Note that the amplitudes of these terms are larger than those detected in the HD 224945 data, giving a higher $`S/N`$ value. The other two terms (1.368, second panel, and 1.145 cd-1, fifth panel) are more separated.

The power spectrum obtained by introducing the five frequencies as k.c.’s is not homogenously flat (bottom panel). The noise in the region 0.0–2.0 cd-1 is 1.34 mmag, while in the region 4.0–6.0 cd-1 it is only 0.65 mmag. However, it should be noted that the general patterns of the residual $`B`$ and $`V`$ power spectra are slightly different. This fact suggests that noise dominates over pulsation and the latter is really undetectable owing to the very small amplitude of the involved terms.

### 5.2 The 1994 and 1991 seasons

Once again, the light curves of the 1991 and 1994 seasons can only be understood if extra alias-free frequency information is available. The frequency analysis of the 1991 $`B`$ data provided evidence for the terms 1.697, 1.627, 1.368 and 1.565 cd-1. The first three terms were already reported in Paper I, where two periodicities having each the shape of a double–wave were proposed. The 1.565 cd-1 (detected also in the 1995 dataset) is a new term obtained by pushing one further step in the analysis. As a result, we also detect here the close triplet observed in the 1995 campaign. Since the frequency resolution in the 1991 dataset is much better than in the 1995 one, this confirmation strengthens our confidence in the reality of the triplet. The peak at 1.145 cd-1 does not appear very clearly in the power spectrum, but it has been considered for the least–squares solution.

The first four terms are also detected in the 1994 $`B`$ and $`V`$ datasets, but residual signal is left, especially in the region 0.5–1.0 cd-1. The 1.145 cd-1 term shows an amplitude much smaller than the other four terms. As a matter of fact, in this season it seems that many terms having amplitudes at the mmag level are excited.

### 5.3 Frequency refinement and least–squares fitting

We followed the same procedure used for HD 224945, refining the frequency values by means of the 1991 dataset and then calculating the amplitudes for the two 1991 subsets and the other datasets. The frequency analysis already gave some hint about the amplitude variability of some terms, since for example the 1.145 cd-1 peak was not clearly detected in all datasets.

This is confirmed by the large amplitude variability, which can be discerned from Tab. 3. Considering the 1991 subsets we can see differences of up to 2.3 mmag, but note the +3.7, +2.9 and –3.9 mmag differences between the amplitudes of the 1.627, 1.368 and 1.565 cd-1 terms in 1995 and 1994 datasets, which look extremely large. The term ranking is also largely changing. Note that the largest amplitude of the 1.565 cd-1 term was in 1994, whereas usually this term was below the first three in 1991 and 1995. The amplitude is also greatly changing. The sum of the squared amplitudes is higher in the 1995 data (311 mmag2) and lower in the 1994 dataset (250 mmag2), with a variation of 25%. This difference supports the intrinsic variability of the mode amplitudes. Note for example that in the 1991 data there are three terms having a large amplitude, in 1994 data there are four and in 1995 data two of them clearly have the largest amplitude, considerably greater than the others. Hence, the amplitude changes look more conspicuous in HD 224638 than in HD 224945.

In the case of HD 224638 we can merge the 1991, 1994 and 1995 datasets knowing that all the datasets show the same frequency content. Indeed, the frequency analysis of the whole set yields the same frequencies listed in Tab. 3. The detection of the three terms at 1.565, 1.627 and 1.697 cd-1 supports the identification as three independent modes proposed in Sect. 5.1. Unfortunately, we cannot give refined values for the frequencies, owing to the cy-1 aliases. However, we performed some tests assuming a constant value for the frequency and then calculated the amplitudes and phases of each term for each observing season. Besides the verification of the amplitude variability, we found that the same frequency displays similar phase values from one season to the next. This fact supports an intrinsic variability of the mode amplitudes, rather than a beating phenomenon between two very close terms.

## 6 Physical properties

The properties of $`\gamma `$ Dor stars have been reviewed by Zerbi (2000); driving mechanisms have been proposed by Guzik et al. (2000), Wu (2002) and Löffler (2002), but the origin of the pulsation in $`\gamma `$ Dor stars still remains an open problem.

HD 224945 and HD 224638 have very similar physical properties (Zerbi 2000): $`M_V`$=2.98, $`L/L_{\mathrm{}}`$=5.5, $`T_{\mathrm{eff}}`$=7200 K, $`R/R_{\mathrm{}}`$=1.51 and $`M/M_{\mathrm{}}`$=1.52 for HD 224638, $`M_V`$=3.07, $`L/L_{\mathrm{}}`$=5.1, $`T_{\mathrm{eff}}`$=7250 K, $`R/R_{\mathrm{}}`$=1.43 and $`M/M_{\mathrm{}}`$=1.51 for HD 224945. Also the difference in metallicity is not significant, \[Me/H\]=–0.15 vs. –0.30, respectively. The two stars occupy very similar positions also in the colour–magnitude diagram, being located on the ZAMS, in the middle of the domain of $`\gamma `$ Dor stars and on the low temperature edge of the $`\delta `$ Sct region (Handler 1999). Therefore, two very similar stars display different pulsational modes, since HD 224945 is characterized by four frequencies spread over a large interval, while HD 224638 displays a much more closely spaced set (it is similar to HR 2740; see Poretti et al. 1997). Moreover, the residual power spectrum of HD 224638 has a very low noise level above 4.0 cd-1 , while that of HD 224945 has a higher level, suggesting some signal contribution.

The large spread in the frequency values observed for HD 224945 is another piece of evidence that the cause of light variability in these stars is pulsation and not stellar activity. Differential rotation is not able to match a spread in frequencies as large as 0.7 cd-1(considering the three terms always observed, i.e., not considering the even more distant frequency at 1.16 cd-1).

We note that multiperiodicity is not observed in all $`\gamma `$ Dor stars. Some display only a single photometric period (i.e., a single low–order mode). A relevant example is HD 207223$``$HR 8330 (Aerts & Kaye 2001), which is also monoperiodic from a spectroscopic point of view (i.e. it does not show any high–order modes). Another example is HD 164615 (Zerbi et al. 1997), one of the first known $`\gamma `$ Dor stars. It is probably a monoperiodic variable which shows amplitude modulation.

Finally, we note that HD 224945 exhibits the shortest known period (i.e. 0.3330 d) of the $`\gamma `$ Dor stars. Guzik et al. (2000) considered 0.48 d as the lower limit (probably the old value of the same term, i.e. 2.00 cd-1, Paper I). Therefore, the period value cannot be used to separate $`p`$– and $`g`$–modes, as $`\delta `$ Sct stars can display such ‘long–period’ $`p`$–modes. A careful evaluation of the physical parameters is necessary. Indeed, in HD 224638 and HD 224945 the fundamental radial mode is shorter than 0.07 d, confirming that $`p`$– and $`g`$–modes are well separated in $`\gamma `$ Dor stars.

## 7 Conclusions

The solution of the light curves of HD 224638 and HD 224945 has been deduced only on the basis of a multisite campaign, since the true peaks could not be recognized if a $`\pm 1`$cd-1 effect is present in the spectral window. We did our best to identify the correct excited modes, but the more important characteristic of these two stars is the strong multiperiodicity itself, independent of the exact frequency values. In our opinion, HD 224638 and HD 224945 provide two of the best examples of features typical of multiperiodicity among $`\gamma `$ Dor variables: different sets of frequency content, amplitude variations, disappearing terms, close doublets of frequencies. Also, after detecting five or six terms, the rms scatter is larger than the observational error. Therefore, residual signal is hidden in the noise. Such extreme multiperiodicity can be clearly ascribed to pulsation, as stellar activity is not able to generate it.

The differences in the frequency and amplitude ranges of HD 224638 and HD 224945 (clearly visible in Fig. 1) remind us of the unpredictable frequency content of $`\delta `$ Sct stars, where the selection mechanism among all the possible modes seems to be different from one star to the next (Poretti 2000). The amplitude variability of the excited modes is another point of similarity between $`\gamma `$ Dor and $`\delta `$ Sct stars. Other than multiperiodicity, HD 224945 and HD 224638 provide the best examples of amplitude variability among $`\gamma `$ Dor stars, following the three campaigns carried out on these stars. This observational evidence suggests strategies for asteroseismological space missions: a long observing run (or two separate runs) may be very helpful in detecting more terms (excited at different levels at different times) or in studying damping effects.

###### Acknowledgements.

ER and SM acknowledge the partial support by the Junta de Andalucia and by the Direccion General de Investigacion (DGI) under project AYA2000-1559. SM also acknowledges the financial support by the Osservatorio Astronomico di Brera and by the Agenzia Spaziale Italiana (ASI Contract I/R/037/01). Thanks are due to an anonymous referee and to M. Breger for useful comments and suggestions on the first version of the manuscript.


text:
# 1 Motivation

## 1 Motivation

Ultra High Energy Cosmic Rays (UHECR) are microscopic particles –protons, photons or perhaps more exotic objects– with a macroscopic energy, about 50 Joules per particle in the extreme of the present observed spectrum. They strike the upper layers of the Earth atmosphere at a rate of about one event per century and per kilometer square. Over one hundred of them have been detected so far by kilometer scale detectors and many more are expected to be seen by forthcoming observatories.

Explaining how these particles get this huge energy is a real challenge for our present understanding of the cosmos. One UHECR primary gets somehow ten million times the energy that a proton will gain at the future Large Hadron Collider (LHC) at CERN. Nature has always been at ease beating the achievements of mankind. One can envisage two broad classes of mechanisms by which nature can produce UHECRs, the so-called bottom up class models and the top-down class models.

In bottom-up models charged particles are accelerated by magnetic fields in large astrophysical sites. In top-down models particles are not accelerated but are created at birth with the huge energy typical of UHECRs.

Here we will concentrate on a particular top-down model . We will assume the UHECRs are produced by the decay of a population of superheavy dark matter particles with a lifetime longer than the age of the Universe $`\tau >10^{10}`$ yr and with mass $`100`$ Joule$`/c^2`$. Theoretical motivation for these particles can be found in .

A superheavy dark matter population created at some stage of the early universe  will gravitationally cluster in the galactic halo. Since the length scale of the halo is around 100 kpc, UHECRs produced in the galactic halo will not have time to interact with the CMB before they reach the Earth. The absence of GZK cut-off is a genuine prediction of models where UHECRs are produced by the decay of a superheavy dark matter particle clustered in the galactic halo.

Since the halo shape is close to spherical one expects a quasi isotropic distribution of events from a halo superheavy dark matter population, which is compatible with experiments. At present it is not possible to make any strong claim about the observed angular distribution of UHECR events because of low statistics. Hopefully, future observatories like Pierre Auger will gather a large enough sample of events to settle down this issue .

Finally there is the question of the UHECR composition. The cosmic ray observatories cannot measure the composition of the primary flux (whether they are photons, protons or heavy nuclei) in an event by event basis. Composition can only be determined in a statistical way. Present analyses tend to favour protons as the main component in the primary flux . In top-down models the main component in the production site, in our case the the galactic halo, is the photon component (actually, neutrinos dominate over photons in a large range of the spectrum but their probability to be detected is too small with present or past detectors). However, on their way to the Earth, photons will interact with the low frequency radio background in the Galaxy and their total flux may be substantially diminished so that on the Earth the total photon flux may be comparable or smaller than the baryon flux. Whether this is indeed possible, subject to the EGRET bound on the low energy $`\gamma `$-rays which will be created by the electromagnetic cascading, will be discussed elsewhere.

Summing up, there are three main tests to falsify or support the hypothesis of UHECRs produced by the decay of dark matter particles clustered in the galactic halo. The first one is the energy distribution of events or spectrum. The second one is the expected angular distribution of events in the sky. The third one is the composition of the primary flux. Here we will focus on the first test and, partially, on the third one. We will briefly show how quantum chromodynamics (QCD) can be used to calculate the spectrum and composition (without photon galactic processing) of the expected UHECR flux. For further details see . QCD was also used to calculate the spectra of UHECRs in .

## 2 DGLAP Evolution

Let us assume that a dark matter particle $`X`$ population, superheavy ($`M_X>3\times 10^{12}`$ GeV) and metastable, with a lifetime larger than the age of the universe, is clustered in the galactic halo. Barring unnatural relations between the coupling constant of $`X`$ and the other fields the total particle multiplicity produced by its decay will be dominated by partonic decays. The flux for the primary $`h`$ (baryon, photon or neutrino) is proportional to the inclusive decay width

$$J^h(E)\frac{1}{\mathrm{\Gamma }_X}\frac{\mathrm{d}\mathrm{\Gamma }(Xh+\mathrm{})}{\mathrm{d}x}=\underset{a}{}_x^1\frac{\mathrm{d}z}{z}\frac{1}{\mathrm{\Gamma }_a}\frac{\mathrm{d}\mathrm{\Gamma }_a(y,\mu ^2,M_X^2)}{\mathrm{d}y}|_{y=x/z}D_a^h(z,\mu ^2).$$
(1)
Particle $`h`$ carries a fraction $`x`$ of the maximum available momentum $`M_X/2`$, and a fraction $`z`$ of the parton $`a`$ momentum. The first term in the integrand is the decay width of $`X`$ into parton $`a`$, $`\mathrm{d}\mathrm{\Gamma }_a/\mathrm{d}y`$, which is calculable in perturbation theory; in lowest order and for 2-body decay it is proportional to $`\delta (1y)`$. The second factor, the non-perturbative $`D_a^h`$, is the fragmentation function (FF) for particles of type $`h`$ from partons of type $`a`$. It gives the expected mean number of particles $`h`$ coming from parton $`a`$. Fragmentation functions cannot be calculated from firsts principles but their dependence in the energy scale $`\mu `$ is governed by the Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) equations .

$$\frac{D_a^h(x,\mu ^2)}{\mathrm{ln}\mu ^2}=\underset{b}{}\frac{\alpha _s(\mu ^2)}{2\pi }P_{ba}(x,\alpha _s(\mu ^2))D_b^h(x,\mu ^2),$$
(2)
where $`\alpha _s(\mu ^2)`$ is the strong coupling constant and $`P_{ba}(x,\alpha _s)`$ is the splitting function for the parton branching $`ab`$. Here the convolution of two functions $`A(x)`$ and $`B(x)`$ is defined as

$$A(x)B(x)_x^1\frac{\mathrm{d}z}{z}A(z)B(\frac{x}{z}).$$
(3)
The splitting functions can be expanded perturbatively:

$$P_{ba}(x,\alpha _s)=P_{ba}(x)+𝒪(\alpha _s).$$
(4)
We limit our study to leading order in $`\alpha _s`$ and therefore ignore $`𝒪(\alpha _s)`$ corrections to the splitting functions. It is also convenient to define the following dimensionless evolution parameter

$$\tau \frac{1}{2\pi b}\mathrm{ln}\frac{\alpha _s(\mu _0^2)}{\alpha _s(\mu ^2)},$$
(5)
$`b`$ being the coefficient in the leading order $`\beta `$-function governing the running of the strong coupling: $`\beta (\alpha _s)=b\alpha _s^2`$. We take $`D_a^h`$ to represent the sum of particle $`h`$ and, if different, its antiparticle $`\overline{h}`$

In the Standard Model there are two partons species: gluons $`g`$ and quarks $`q`$ (for the sake of brevity we shall only consider the singlet quark, i.e., the sum of all quark and antiquark flavours). The DGLAP equations can be written as

$$_\tau \left(\begin{array}{c}D_q\hfill \\ D_g\hfill \end{array}\right)=\left(\begin{array}{cc}P_{qq}& 2n_\mathrm{F}P_{gq}\\ P_{qg}& P_{gg}\end{array}\right)\left(\begin{array}{c}D_q\hfill \\ D_g\hfill \end{array}\right).$$
(6)
When supersymmetry (SUSY) is included one has in addition the gluinos and the squark singlet. The DGLAP equations are now

$$_\tau \left(\begin{array}{c}D_q\hfill \\ D_g\hfill \\ D_s\hfill \\ D_\lambda \hfill \end{array}\right)=\left(\begin{array}{cccc}P_{qq}& 2n_\mathrm{F}P_{gq}& P_{sq}& 2n_\mathrm{F}P_{\lambda q}\\ P_{qg}& P_{gg}& P_{sg}& P_{\lambda g}\\ P_{qs}& 2n_\mathrm{F}P_{gs}& P_{ss}& 2n_\mathrm{F}P_{\lambda s}\\ P_{q\lambda }& P_{g\lambda }& P_{s\lambda }& P_{\lambda \lambda }\end{array}\right)\left(\begin{array}{c}D_q\hfill \\ D_g\hfill \\ D_s\hfill \\ D_\lambda \hfill \end{array}\right).$$
(7)
We have developed a C++ code to solve the DGLAP equations for fragmentation functions. We have extended the algorithm introduced in Ref.  to be able to solve the SUSY equations as well as the Standard Model equations .

We begin with FFs at the energy scale $`\mu _0=M_Z`$ and evolve them to the final energy scale $`\mu =M_X`$ using Eq. (6) and also Eq (7) if SUSY is included. For the initial fragmentation function of baryons, $`D_a^p(x,M_Z^2)+D_a^n(x,M_Z^2)`$, we adopt the fit performed in Ref.  to LEP hadronic data . For photons and neutrinos we generate initial data at the $`Z`$ peak using the QCD Monte Carlo event generator HERWIG . Comparison with LEP data shows that although HERWIG overproduces baryons at high $`x`$ , its photon and meson output at the $`Z`$ peak matches the experimental spectra remarkably well. Since neutrinos mainly come from charged pion and kaon decays, one can thus be confident in taking the HERWIG generated FF for neutrinos as the initial condition for the evolution. There is also a sizable contribution from heavy flavour decays to the neutrino spectrum at high $`x`$ which is explicitly taken into account by HERWIG.

A (s)parton is not included in the evolution as long as the energy scale is lower than its mass; when the threshold for its production is crossed, it is added to the evolution equations with an initially vanishing FF and it is assumed to be a relativistic particle.

In the SM case we evolve the $`q`$ and $`g`$ initial fragmentation functions from $`M_Z`$ to $`M_t`$, the top quark mass, with the number of flavours set to $`n_\mathrm{F}=5`$, and then evolve from $`M_t`$ to $`M_X`$ with $`n_\mathrm{F}=6`$.

In the SUSY case we evolve the $`q`$ and $`g`$ initial fragmentation functions from $`M_Z`$ to the supersymmetry breaking scale $`M_{\mathrm{SUSY}}>M_t`$ using the SM equations to obtain $`D_i^h(x,M_{\mathrm{SUSY}}^2)`$, with $`i=q,g`$. Then we take $`D_i^h(x,M_{\mathrm{SUSY}}^2)`$, $`i=q,g`$, and $`D_j^h(x,M_{\mathrm{SUSY}}^2)=0`$, $`j=s,\lambda `$, and evolve them from $`M_{\mathrm{SUSY}}`$ to $`M_X`$ using the SUSY equations. All spartons are taken to be degenerate with a common mass $`M_{\mathrm{SUSY}}`$.

## 3 Ultra High Energy Cosmic Ray Spectrum

We can now translate the calculated fragmentation functions into the expected cosmic ray spectrum in order to confront the observational data. In the previous Section (2) we have briefly shown how to calculate quark singlet and gluon functions for the SM, and quark singlet, gluon, squark singlet and gluino functions for SUSY. In the absence of a specific model for the different branching ratios we weight all the (s)parton contributions evenly. For 2-body decay of $`X`$, $`x=2E/M_X`$ the flux of particle $`h`$ is given by :

$$E^3J^{\mathrm{halo}}(E)=Bx^3D^h(x,M_X^2).$$
(8)
We have multiplied the flux by $`E^3`$, as is usual, to emphasise the structure in the spectrum near the GZK energy. The normalisation factor $`B`$ is common for the galactic halo flux of baryons, neutrinos and photons, and determines the quantity $`n_X/\tau _X`$ (see Eq. 1).

Let us now compare the calculated cosmic ray flux to the published data from Fly’s Eye , AGASA , Haverah Park  and Yakutsk . In Ref.  these data have been carefully assessed for mutual consistency and appropriate adjustments made to the energy calibration. Its authors recommend adoption of the following standard differential energy spectrum below the GZK energy, in the range $`4\times 10^{17}\mathrm{eV}<E<6.3\times 10^{18}\mathrm{eV}`$:

$$J(E)=(9.23\pm 0.65)\times 10^{33}\text{m}^2\text{s}^1\text{sr}^1\text{eV}^1\left(\frac{E}{6.3\times 10^{18}\text{eV}}\right)^{3.20\pm 0.05},$$
(9)
with the spectrum flattening at higher energies as $`J(E)E^{2.75\pm 0.2}`$ up to the GZK energy, and extending further to at least $`3\times 10^{20}\mathrm{eV}`$ . Thus the UHECR spectrum can naturally be interpreted  as the superposition of the ‘low energy’ component (9), and the new ‘flat’ component that extends into the post-GZK region. The former is presumably galactic in origin (consistent with the detection of anisotropy at $`10^{18}\mathrm{eV}`$ ), while the latter is interpreted  as produced by the decay of a superheavy particle population in the galactic halo. Taking baryons to be the dominant primary UHECRs as indicated by experiment, the total flux is

$$E^3J(E)=\frac{k}{E^m}+Bx^3D^{\mathrm{baryon}}(x,M_X^2),$$
(10)
where the values of $`k`$ and $`m`$ can be read off Eq. (9). Note that since $`D^{\mathrm{baryon}}`$ and $`D^\gamma `$ have a similar shape, taking photons to be the primaries would just alter the normalisation $`B`$. In Fig. 1 we plot the best SM evolution fit to the cosmic ray data while in Fig. 2 we plot the best SUSY evolution fit.

The assumption of 2-body decay may be rather naive for a superheavy particle like a crypton  which is expected to decay through very high-order non-renormalisable operators. Many-body decay distributes the total energy $`M_X`$ among several particles and thus flattens the spectrum. We assume that many-body effects are purely kinematical and hence can be encapsulated in the phase space of the decay. Let $`\rho _n(z)`$ be the probability density that one parton carries off a fraction $`z`$ of the total available energy per parton $`M_X/2`$. For $`n3`$ we get

$$\rho _n(z)z(1z)^{n3}.$$
(11)
To leading order in $`\alpha _s`$ the particle flux is then given by

$$E^3J^{\mathrm{halo}}(E)=Bx^3_x^1\frac{\mathrm{d}z}{z}\rho _n\left(\frac{x}{z}\right)D^h(z,M_X^2).$$
(12)
In particular if $`D^h(x,M_X^2)(1x)^{a(M_X^2)}`$ as $`x1`$, the differential particle flux decreases as $`J^{\mathrm{halo}}(x)(1x)^{a(M_X^2)+n2}`$.

For many-body $`X`$ decays the total flux is

$$E^3J(E)=\frac{k}{E^m}+Bx^3_x^1\frac{\mathrm{d}z}{z}\rho _n\left(\frac{x}{z}\right)D^{\mathrm{baryon}}(z,M_X^2),$$
(13)
where $`n`$ is the number of partons into which $`X`$ decays. In Fig. 3 we plot evolved SUSY spectra for different values of $`n`$.

## 4 Conclusions

We have calculated the UHECR spectra of baryons, photons and neutrinos expected from the decay of superheavy dark matter particles of mass $`M_X10^{12}`$ GeV which are clustered in the galactic halo. The spectrum for every primary particle is given by FFs at the scale $`M_X`$. We have calculated these FFs using QCD DGLAP equations. We have considered Standard Model equations and SUSY equations. We have also taken into account the possibility of many-body decay in the decay of $`X`$. The shape of the fragmentation spectrum (of either baryons or photons) fits rather well the new component of ultra-high energy cosmic rays extending beyond the GZK energy.

ACKNOWLEDGEMENTS

The work here outlined was performed in collaboration with Subir Sarkar. I would like to acknowledge support by a Marie Curie Fellowship No. HPMF-CT-1999-00268.