© ESO 2018

1 Introduction

Young stars of spectral types B and O are frequently found in binary and multiple systems (for a review, see Duchêne & Kraus 2013). Recently, it has become clear that the preference for very close binaries is stronger than generally assumed. Especially for stars at the high-mass end of the spectrum, interaction with a binary companion appears to be the dominant mode of evolution (Kobulnicky & Fryer 2007; Kiminki & Kobulnicky 2012; Sana et al. 2012; Chini et al. 2012; Almeida et al. 2017). In galaxies, these O-type and B-type stars and their more evolved counterparts play crucial roles in chemical evolution and feedback. They are the progenitors of white dwarfs, neutron stars, and black holes that can give rise to many exotic phenomena, for example explosive transients, especially when they interact with a binary companion. Improving our understanding of the products of binary interaction is therefore relevant for multiple disciplines in astrophysics including the new field of gravitational wave astrophysics.

A large fraction of close binary systems evolve through a relatively simple first post-interaction phase. This phase starts after the (first) episode of stable Roche-lobe overflow ceases. The original donor star has then lost nearly its entire hydrogen-rich envelope. The star becomes hot and compact and in most cases resides close to the helium main sequence (e.g., Kippenhahn & Weigert 1967; Pols et al. 1991; Yoon et al. 2010). The less evolved companion still resides on the main sequence. Typically it now is the brightest star in the system, at least in the optical part of the spectrum. The stripped star emits most of its light in the far- or extreme-UV (Götberg et al. 2017). This phase is usually relatively long-lived as both stars evolve on the slow nuclear burning timescale. We refer to systems in this post-interaction phase as stripped-envelope star (SES) systems. The only way to avoid this phase is if the stars merge (e.g., Podsiadlowski et al. 1992) or experience certain proposed minority channels, for example those involving chemically homogeneous evolution (de Mink et al. 2009), as proposed as a formation channel for binary black holes (Mandel & de Mink 2016; Marchant et al. 2016). The SES phase may be very short-lived for very wide systems where the donor fills its Roche lobe after completion of central helium burning (case C mass transfer).

Contrary to what might be expected for a common long-lived phase, we only have sporadic identifications of SES systems. This poses an apparent paradox: despite the high binary fraction among young stars, these post-interaction binary systems containing an SES are rarely observed. This may simply be the result of observational biases, as argued for example by de Mink et al. (2014) and Götberg et al. (2017). Post-interaction systems are very hard to detect due to a variety of selection effects and mostly pose as apparently single stars. While this may provide a solution to the paradox, the situation remains unsettling. Further insight into post-interaction systems and their detectability is warranted, especially because they can provide stringent tests for the uncertain mass transfer process.

The few observed systems that contain a main sequence star + SES (e.g., the Be + sdO star systems; see Sect. 6) clearly show how hard it is to detect them. The high temperatures of SESs hinder their detection due to the lack of absorption lines resulting from a high degree of ionization. Also, because of their high temperatures SESs shine most brightly in the UV; the companion typically outshines them in optical bands, which are more accessible to observations (e.g., Götberg et al. 2017). Detecting their excess at short wavelengths requires observations from space, for example with the Hubble Space Telescope. These facilities arescarce and highly oversubscribed. Indirect evidence for the presence of these excesses can be obtained if the main sequence companion shows radial velocity variations due to the orbital motion. Unfortunately, these variations are often small because of an extreme mass ratio after mass transfer. Moreover, the main sequence star is expected to rotate rapidly (e.g., de Mink et al. 2013). Line broadening due to rotation and possible variations in the line profile due to the presence of an outflowing disk complicate the detection of small radial velocity variations.

The Be + sdO binary system φ Persei is one of the few exceptions where the stripped remnant of a donor star has been detected. We elaborate on the observational history of this system in Sect. 2. The sdO star is thought to be stripped of its envelope by the now rapidly rotating Be star (e.g., Gies et al. 1998), as has also been shown based on dedicated theoretical simulations for this system by Vanbeveren et al. (1998a,b) and Pols (2007). Recently, new data for this system have been presented by Mourard et al. (2015) who unambiguously verified the presence of the SES by directly imaging the sdO component of φ Persei using new high angular resolution observations. Their astrometric data led to a slightly revised orbital solution, but are consistent with those derived earlier by Gies et al. (1998). The new observational data motivated us to reinvestigate the evolutionary status of the φ Persei system and expand upon the theoretical studies presented earlier by Vanbeveren et al. (1998a,b), and Pols (2007) to find constraints on the physical processes. In particular, we investigate the surprisingly high luminosity of the stripped star (its main sequence companion is eight times more massive, but only two times more luminous) and discuss the implications for its evolutionarystatus. We also place it in the context of analog systems and use it to derive clues to better understand the apparent paradox of the scarcity of SESs.

The structure of this paper is as follows. In Sect. 2, we describe the observational history and characteristics of the φ Persei system. To investigate the φ Persei system we compute a large grid of binary evolutionary models with the fast binary evolution code binary_c and selected models with the detailed evolutionary code MESA, as is described in Sect. 3. In Sect. 4 we present our main results. In Sect. 5 we discuss model variations and alternative implications, and compare them with previous studies. In Sect. 6 we give an overview of φ Persei analogs that have been observed in different evolutionary phases. Section 7 summarizes our conclusions.

2 Observations of the φ Persei system

At a distance of 183 ± 3 pc and an apparent magnitude of m  =  4.09 (Mourard et al. 2015, and references therein), φ Persei can be observed with the naked eye. The main sequence star itself is classified as Be star of type B1.5 V:e-shell (Slettebak 1982): the “e” indicates spectral emission features and “shell” indicates that both rotationally broadened lines and narrow absorption lines are present, the second presumably originating from a circumstellar shell. At the beginning of the 20th century, radial velocity variations were already being detected and directly interpreted as a sign of binarity (Campbell 1902; Ludendorff 1910; Cannon 1910).

The nature of the secondary has been the matter of long-standing debate. Poeckert (1979, 1981) reported the discovery of a weak He II λ4686 Ångström emission line and found that it showed radial velocity variations moving in antiphase with the lines of the Be star primary. He argued that the emission line originated from hot gas surrounding the unseen companion. To explain the hard spectrum of the companion needed to ionize helium and produce the helium emission line, the author proposed that the companion could be a remnant core of a star stripped in a mass transfer event. Gies et al. (1993) discovered emission components of the He I λ6678 Ångström line and found that they also moved in antiphase with the lines of the primary, supporting the findings by Poeckert (1981). An overview of the available data at that time and a revised orbital solution was presented by Bozic et al. (1995).

Motivated by the prediction that the still unseen companion would be bright in the UV, Thaller et al. (1995) used the International Ultraviolet Explorer to search for direct features of the companion. They found that the spectrum of the primary shows evidence of rotationally broadened photospheric lines, consistent with the expectation for a rapidly rotating Be star. They reported the discovery of many weak photospheric lines similar to those found in hot subdwarf O (sdO) stars, as well as strong emission of C IV λ1550 Ångström. Even though the spectrum was noisy (the hot companion only contributes about 12% of the UV flux), it provided the first direct confirmation of the nature of the companion. Gies et al. (1998) were able to confirm their findings by clearly showing the features of the faint hot subdwarf companion, using the Goddard High Resolution Spectrograph, which was one of the original instruments on board the Hubble Space Telescope.

Recently, Mourard et al. (2015) obtained the full orbital solution of the φ Persei system through high angular resolution imaging observations. To achieve this, they used long-baseline interferometry in the near-IR and optical, and then combined their astrometric data with the earlier obtained radial velocity measurements. According to the latest knowledge, the φ Persei system consists of a very rapidly rotating (v_rot = 0.93 ± 0.08 v_crit) Be star withmass M = 9.6 ± 0.3M_⊙ and a hot SES of M = 1.2 ± 0.2M_⊙ with a 126-day orbital period, consistent with the parameters obtained earlier by Gies et al. (1998). Despite its much lower mass, the SES is almost half as luminous (L_SES = 10^3.8±0.13L_⊙) as its Be star companion (L_Be = 10^4.16±0.10L_⊙). The derived parameters of the φ Persei system are listed in Table 1.

3 Method

To explore the parameter space with a large grid of evolutionary models we use the rapid synthetic binary evolution code binary_c (Sect. 3.1). The advantage of this code is that it is fast enough to synthesize a large population of stars to use for our analysis. The downside is that certain physical processes are treated in less detail than in a full stellar evolutionary code. Therefore, we use the detailed stellar evolution code MESA (Sect. 3.2) to test the evolutionary tracks of SESs that are predicted with binary_c.

3.1 Synthetic binary evolution with binary_c

This code was originally developed as a rapid single-star evolution code (Hurley et al. 2000) based on analytical fits to the detailed stellar evolution models provided by Pols et al. (1998). Later, it was made suitable to follow the evolution of binary systems by including physics such as tidal interaction and mass transfer, as described in Hurley et al. (2002). We use a version of the code known as binary_c developed and updated by Izzard et al. (2004, 2006, 2009) and de Mink et al. (2013) to include more up-to-date physics assumptions. We adopt the same assumptions as listed in de Mink et al. (2013, and references therein), except that we ignore rotationally enhanced mass loss via stellar winds. A full description can be found in thereferences above; here we give a summary of the most relevant physics assumptions that are in this work.

Overshooting is accounted for in the detailed evolutionary models that binary_c relies on, using a parametrization for the overshooting parameter as described in Pols et al. (1998). In our mass range of interest it corresponds roughly to a classical overshooting parameter α_OV ≈ 0.28, which refers to the extent of the mixed region above the convective core expressed in pressure scale heights. The assumed overshooting parameter is relevant since it determines the size of the convective core during the main sequence evolution, and thus the mass of the resulting helium core. This in turn sets the mass of the stripped star after mass transfer.

Mass transfer and the criteria for stability of Roche-lobe overflow are modeled as described in Hurley et al. (2002). For post-main sequence donors, mass transfer is assumed to lead to complete stripping of the hydrogen envelope. The further evolution is computed using fits against evolutionary tracks for naked helium star models (Pols et al. 1998).

We account for the rejuvenating effect of mass accretion using the algorithms described in Hurley et al. (2002) but using an updated treatment for the core masses, as described in de Mink et al. (2013). This reduces the apparent age of the star while its mass increases, resulting in an increase in the temperature and luminosity of the accretor. We also follow the angular momentum evolution of both stars and model the effect of spin up of the mass gainer as described in de Mink et al. (2013). We do not account for the possible additional effect of rotationally induced mixing processes or gravity darkening. We discuss possible concerns resulting from limitations in the treatment of rotation in Sect. 4. These affect the comparison with the rapidly rotating Be star. Our predictions for the stripped star, the main focus of this study, are not significantly affected by these limitations.

To allow for non-conservative mass transfer, we limit the accretion rate Ṁ_acc to ten times the mass of the accretor star M_a divided by its Kelvin–Helmholtz timescale τ_KH (the time a star needs to radiate away its potential energy). If this limit is not exceeded, we assume that all mass lost from the donor star is accreted by its companion. Therefore we have $${\dot{M}}_{\text{acc}}=\text{min}\left(-{\dot{M}}_{\text{donor}},\frac{10M_{\text{a}}}{{\tau }_{\text{KH}}}\right).$$(1)

To investigate the impact of uncertainties in the mass transfer efficiency β = Ṁ_acc∕Ṁ_donor, we also consider cases where we vary β ∈{0.25, 0.5, 0.75, 1.0} following de Mink et al. (2007). We thus implicitly assume that the accretion disk regulates the angular momentum that is accreted by the mass gainer (see discussion and Appendix A.3 in de Mink et al. 2013). We assume that all mass that is not accreted is lost from the system taking away the specific angular momentum of the orbit of the accretor star (e.g., Petrovic et al. 2005; van den Heuvel et al. 2017).

3.1.1 Synthesized population

With binary_c, we evolve a grid of binaries until both stars in each system have left the main sequence. The following prior (“birth”) distribution functions are assumed:

• For the primary mass distribution: $$f_{{\text{M}}_1}(M_1)d{M}_1\propto M_1^{\alpha }d{M}_1.$$(2)

  The chosen mass range is 4.5M_⊙ ≤ M₁ ≤ 9.0M_⊙, which is large enough to find a solution for the φ Persei system. An initial M₁ lower than ~5.5M_⊙ is ruled out because the total mass of the system is ~11M_⊙ and M₂ < M₁ by definition. The upper limit is based on a rough grid with a larger parameter space that was synthesized earlier, which showed no solutions at higher masses. For this mass range we assume α = −2.35, ergo a Salpeter initial mass function (Salpeter 1955).

• The mass of the secondary star is drawn from the mass ratio (q = M₂∕M₁) distribution: $$f_q(q)dq\propto q^{\kappa }dq.$$(3)

  This distribution is generally believed to be flat (see, e.g., Kobulnicky & Fryer 2007; Kouwenhoven et al. 2007). Therefore, we use κ = 0. We explore the mass ratio range 0.1 ≤ q ≤ 1.0.

• The chosen period distribution is flat in logspace, also known as Öpik’s law (Öpik 1924), $$f_P(P)dP\propto \frac{1}{P}dP,$$(4)

  since 1∕P dP = dlogP. The chosen period rangeis 2.5 days ≤ P ≤ 40 days. The limits are again based on the rough grid with a larger parameter space, which showed no solutions at lower and higher periods.

A flat prior probability distribution for the age of the system is used for the analysis. The size of the grid is 24 × 24 × 24 (N(M₁) × N(q)  × N(P), 13824 systems in total). We only consider circular systems; Mourard et al. (2015) find that the current orbital solution converges a system with zero eccentricity. We assume a metallicity of Z = 0.02 for the population. This metallicity matches the derived metallicity of the α Persei cluster (Netopil & Paunzen 2013) with which φ Persei is associated (see Sect. 4.2).

3.1.2 Bayesian analysis

We adopt a Bayesian approach that uses the prior probability of initial parameters to extract the posterior probability distribution functions using the synthesized population. Bayes’ theorem can be expressed as $$P(M\mid D)=\frac{P(D\mid M)P(M)}{P(D)},$$(5)

where M refers to a model parameter and D refers to the observed data, and thus P(M∣D) is the posterior probability distribution of an M given the data from the observations (D). Here, we are interested in obtaining the posterior probability distributions of the following model parameters M: the birth masses of both stars M_SES,0 and M_Be,0, the initial orbital period P_orb,0, and the system age. P(M) is a model parameter’s prior probability function, e.g., the initial mass function (Eq. (2)). P(D) is not dependent on model parameters and serves as a normalization factor.

The likelihood distribution of data points for a given a model parameter is expressed as P(D∣M). The observed parameters that we use to determine the likelihood of fit P(D∣M) of the model parameters are the masses and the temperatures of both stars in the binary, the orbital period, and the ratio of the radii, which is a well-constrained value for the system. Therefore, we use six fit parameters in total (M_SES, M_Be, T_{eff, SES}, T_{eff, Be}, P_orb, and R_SES∕R_Be) for the case of φ Persei. Their valuesand errors are given in Table 1. The error on the observed orbital period is extremely small with σ_P,obs = 0.0071 days in comparison to the rather coarse spacing in the initial orbital period in our model grid. To prevent the orbital period from dominating the solution, we adopt instead an error of 15 days. This corresponds to a relative error of about 10%, which is similar to the relative errors on the other parameters (cf. de Mink et al. 2007).

In practice, the method used to obtain P(M∣D) is the following:

• The likelihood of the fit to the observed parameters is calculated for each time step in each simulated binary system using theχ² method. This corresponds to P(D∣M), and it is obtained using $${\chi }^2={\displaystyle \sum _i{\left(\frac{x_{\text{i},\text{obs}}-x_{\text{i},\text{mod}}}{{\sigma }_{\text{i},\text{obs}}}\right)}^2},$$(6)

  where x_i,obs and σ_i,obs are the observed value and error for the fit parameter, respectively. The value of the fit parameter predicted by the model is represented by x_i,mod.

• The likelihood of the fit in each time step is then multiplied by the prior probability of the system, P(M), which results from the probability distributions of the initial parameters M_1,0, q₀ and P_orb,0, which are provided in Sect. 3.1.1.

• The value obtained at step 2 is multiplied by the time step and binned for the initial parameters M_1,0, M_2,0, and P_orb,0 and for the system age.

The obtained result represents the posterior probability distribution of the system’s initial parameters P(M∣D). In other words this is the probability distribution of the birth parameters. We note that the error bars that we derive propagate from the observational constraints and do not take into account model uncertainties.

3.2 Detailed binary evolution with MESA

The main focus of this study is the SES (i.e., the sdO star) in the φ Persei system. Therapid synthetic evolutionary code described above treats these stars in a simplified way as pure helium stars. Using detailed evolutionary simulations (Götberg et al. 2017) has shown that these stripped stars typically still have a thin remaining layer of hydrogen (see also Yoon et al. 2017). To investigate the consequences of this we also use the detailed binary evolutionary code MESA (version 7624, Paxton et al. 2011, 2013, 2015) to produce more realistic models for SESs. This code solves for the one-dimensional stellar structure equations and composition changes simultaneously, and follows the evolution of binary systems as they evolve through mass transfer.

We adopt the Ledoux criterion for convection, while using semiconvection with a value for the efficiency parameter of α_SC= 1 (Langer et al. 1983). Mixing in the convective regions is described by mixing length theory (Böhm-Vitense 1958) and for the mixing length parameter we use α_MLT = 2. We adopt an overshooting parameter α_OV = 0.28 to be consistent with the models by Pols et al. (1998) that are used in our rapid synthetic code binary_c. We also take into account the effects of thermohaline mixing and rotation. The metallicity is again set to Z = 0.02. We use the nuclear network approx21, which contains 21 isotopes and is appropriate for modeling past central helium exhaustion.

We account for mass loss through stellar winds, although the effect is very small in the mass range of interest. We use the algorithm by Vink et al. (2001) for stars with T_eff> 10 kK and a surface hydrogen fraction of X_H,s > 0.4. For hot stars that have lost part or most of their hydrogen-rich envelope, T_eff> 10 kK and X_H,s < 0.4, we adopt the algorithm by Nugis & Lamers (2000). For cool stars, T_eff< 10 kK, the prescription from de Jager et al. (1988) is used. During stable Roche-lobe overflow we use the implicit Ritter scheme (Ritter 1988) to compute the mass loss rate from the donor star. Mass transfer is non-conservative, as described in Götberg et al. (2017), where mass loss from the system is assumed to take away the angular momentum of the accreting star. We do not use MESA to further investigate the evolution of the accretor.

We simulate a binary system which has the same initial orbital period (16 days) and initial secondary mass (3.8 M_⊙) as the best-fitting system we derive with binary_c. The initial primary mass is chosen such that the SES has the desired mass after envelope stripping.

4 Results

Here, we discuss the results of our Bayesian analysis of the φ Persei system against our grid of binary evolutionary tracks. To test our solution, we first describe the posterior distributions of fit parameters and compare them with the observed parameters (Sect. 4.1) followed by a discussion of the properties of the progenitor system φ Persei (Sect. 4.2) and its inferred age (Sect. 4.3). We then discuss the current evolutionary stage (Sect. 4.4) and its future evolution predicted by our best fit model (Sect. 4.5).

4.1 Posterior distributions of the observed fit parameters

In this section, we consider the posterior probability distribution of the fit parameters to test how well the models reproduce the observed system parameters. In an idealized case, both would be identical. Figure 1 shows a comparison of the observed parameters and their 1 σ uncertainty intervals with the posterior probability distribution of the fit parameters obtained from our analysis. The posterior probability distributions are nicely represented by a simple Gaussian, which indicates that they converge to a unique solution.

Generally, we find good agreement between the observed parameters and the posterior probability distributions, which all peak within 1 σ of the observed value. The only exception is the effective temperature of the Be star, T_{eff, Be}. Our best-fitting models have temperatures that are 1.5 σ colder than the observed temperature of the Be star. The difference may be of physical origin. It may indicate that the Be star has experienced extra mixing as a result of mass transfer. Stars that have a more homogeneous interior composition profile tend to be more compact and hotter. Such mixing could for example result from thermohaline mixing if the Be star accreted helium rich material (Stancliffe & Glebbeek 2008) or rotational mixing (Cantiello et al. 2007) as a result of spin-up. In this case, these effects would more than compensate for gravity darkening (von Zeipel 1924). We do not account for these processes in our simulations; however, this does not significantly affect our findings for the evolutionary phase of the stripped object (see Sect. 4.4), which is the main focus of this work.

Fig. 1 Posterior probability distributions of the fit parameters (yellow histograms) of the φ Persei system together with a fitted Gaussian curve (in black). For comparison, the observed parameters for the φ Persei system (bluevertical dashed lines) and their 1 σ uncertaintyintervals (blue shaded regions) are shown. Top panels: current masses of the sdO and Be star, MSES and MBe. Central panels: effective temperatures, Teff, SES and Teff, Be. Bottom panels: present-day orbital period Porb and the ratio of their radii RSES∕RBe.

4.2 Constraints on the progenitor system of φ Persei

In Fig. 2, we show the posterior probability distribution functions for the birth parameters derived using our default simulations. These birth parameters are the initial masses of both stars and the initial orbital period. We also show the current age of the system. The binned distributions are again nicely described by a simple Gaussian fit.

We find that the most probable initial masses of the system are 7.2 ± 0.4M_⊙ for the primary star and 3.8 ± 0.4M_⊙ for the secondary star. Interestingly, the present-day mass of the Be star is just above the threshold for a single star to end its life in a core-collapse supernova (~8M_⊙), while its birth mass was below this threshold. This means that accretion truly changed the final fate of the secondary star in this system, which would have been destined to end its life as a white dwarf if it had been single.

For this initial orbital period, 16 ± 4 days, the system was wide enough to remain detached during the main sequence evolution of the primary star. Mass transfer started after the initial primary started crossing the Hertzsprung gap, which is called case B mass transfer.

We note that our findings are slightly different from but consistent with those derived earlier by Vanbeverenet al. (1998a) and Pols (2007). They both inferred initial masses for the progenitors of around 6M_⊙ for the SES and 5M_⊙ for the Be star, and an initial orbital period on the order of 10 days. The differences are due in part to the fact that they tried to fit an sdO star mass of 1.14 ± 0.04M_⊙ (Gies et al. 1998), which was the best estimate available at the time, while we use the more recently determined value of 1.2 ± 0.2M_⊙ instead. A further difference is that we simultaneously fit for multiple observed parameters at once. We do not fit the luminosity directly since this is a dependent parameter, but by fitting for the temperatures and ratio of the radii a fit solution that can explain the high luminosity of the sdO star is also favored, a solution that prefers a higher mass.

Fig. 2 Posterior probability distributions of the birth parameters (green histograms) of the φ Persei system together with a fitted Gaussian curve (in black). Top panels: birth masses of the sdO star and the Be star, MSES, 0 and MBe, 0. Bottom panels: initial orbital period Porb, 0 and the system age.

4.3 Age of φ Persei: a blue straggler of the α Persei cluster?

The age we derive for the φ Persei binary system in the standard simulation is 57 ± 9 Myr (Fig. 2). Mourard et al. (2015) argue that φ Persei is likely a member of the α Persei cluster because its distance, kinematics, and angular position are comparable to those of the cluster. Our age estimate is in excellent agreement with recent estimates of the age of the α Persei cluster, which is in the vicinity of φ Persei: these range from 52 Myr (Makarov 2006) to 60 Myr (Zuckerman et al. 2012).

This age corresponds to a turn-off mass of about ~ 7M_⊙ star, which we consider the most likely initial mass of the original primary star in the φ Persei system. This by itself supports the picture that the φ Persei Be star is a member of the cluster and that it gained a significant amount of mass from its companion later in its evolution and then rejuvenated. The Be star can be considered a rejuvenated blue straggler member of the cluster.

We note that the cluster age is debated. For example, Stauffer et al. (1999) found a cluster age of 90 ± 10 Myr using a method based on lithium depletion in brown dwarfs. This is not consistent with our findings. It is also on the high side of the highest lifetime that the φ Persei system could have given the minimum mass of the primary star (5.5M_⊙, i.e., half the total system mass), which has a main sequence lifetime of ~80 Myr.

Fig. 3 Hertzsprung–Russell diagram showing the evolutionary track of a stripped envelope star from the best-fitting system derived with binary_c. This object has a mass of 1.36M⊙. It is simulated with both binary_c (dashed line) and MESA (solid line). Models in the helium core burning phase are displayed in blue, models in the helium shell burning phase in orange.

4.4 Evolutionary stage of φ Persei: caught in the short-lived helium shell burning phase

We find that the best fit models for φ Persei imply that the stripped star resides in the helium shell burning phase. This is illustrated in Fig. 3, where we compare thepresent-day parameters of the sdO star with the evolutionary track of our best fit model for the stripped star. We only show the stripped star phase after the first Roche-lobe overflow has ceased. This includes the phase during which the star undergoes central helium burning. After exhaustion of helium in the core, burning continues in a shell. During this phase the star becomes more luminous and swells up.

In our binary_c simulation the stripped star is treated as a pure helium star. For comparison we also show the corresponding evolutionary track of a stripped star resulting from an MESA simulation. The MESA model is slightly cooler during the helium core burning phase. This is due to a thin layer of hydrogen that is left at the surface, which sustains very weak burning of hydrogen in a shell around the core (Götberg et al. 2017). The binary_c simulations ignore this thin hydrogen layer. However, the differences between the two tracks are small; they are comparable to the quoted observational error bars. We therefore conclude that our rapid binary_c simulations are adequate for our analyses.

Figure 3 thus indicates that the best-fitting models for the φ Persei SES assumethat it resides in the helium shell burning phase. The present-day luminosity is nearly 1 dex brighter than the luminosity expected for models in the the helium core burning phase. This is true not only for the best fitting model, but it also holds for the full ensemble of likely solutions for φ Persei. We find a posterior chance of 99.9% that the system resides in the shell burning phase.

Figure 4 further illustrates this by comparing the observationally derived present-day mass and luminosity with the range of luminosities expected for helium core burning stars and helium shell burning stars as a function of mass. For reference we show evolutionary tracks of 0.8, 1.2, and 1.6M_⊙ SESs simulated with both MESA and binary_c. The effect of stellar wind mass loss is visible for the MESA, but apart from this we find that the differences are minimal. For both sets of models it is evident that φ Persei is almost an order of magnitude too luminous to be helium core burning. Even if the mass of the object is 2 σ higher than the observationally derived value, the models strongly favor a helium shell burning solution.

To put the apparently high luminosity of the φ Persei SES into perspective, we also show the parameters for the two other known sdO + Be binaries, FY CMa (Peters et al. 2008) and 59 Cyg (Peters et al. 2013) in Fig. 4. Their observed parameters are summarized in Table 3.

Although it has the same mass within error bars (1.3 ± 0.2M_⊙) as φ Persei (1.2 ± 0.2M_⊙), the FY CMa SES is about 0.5 dex less luminous (Fig. 4), which by itself is already surprising. Because of this it fits to thehelium core burning models and to the helium shell burning models (Fig. 4). It is remarkable that this object has been detected at all given that its luminosity is an order of magnitude lower than the luminosity of its Be star companion (Table 3). A possible explanation may lie in the extensive amount of UV observation time that has been dedicated to this system: between 1979 and 1995, 97 spectra were obtained with the International Ultraviolet Explorer satellite (see Peters et al. 2008, and references therein).

The SES in 59 Cyg has a lower mass (0.77 ± 0.15M_⊙) than the other two SESs. It seems, similar to the φ Persei SES, to be too luminous for its mass to be in the helium core burning phase by about 0.6 dex, which is over 2 σ (Fig. 4). This system seems to be a lower mass counterpart of φ Persei, both being in the shell burning phase.

In summary, we conclude that the φ Persei SES is remarkably bright for its mass. It is about 0.5 dex brighter than the stripped star in FY CMa, which has the same mass within error bars. This is something we can only reconcile in our models if it resides in the short-lived helium shell burning phase.

Fig. 4 Evolution of the luminosity and mass of stripped envelope stars of 0.8M⊙, 1.2M⊙, and 1.6M⊙. Helium core burning models are displayed in blue, while helium shell burning models are displayed in orange. The blue shaded area indicates the parameter space in which the models predict SESs to be helium core burning; in the area above, SESs are predicted to be helium shell burning. In the gray hatched region, the models do not predict helium burning SESs. Top (solid lines): MESA tracks. Bottom (dashed lines): binary_c tracks.

4.5 The future of the φ Persei system: a possible progenitor of a late core collapse or exotic merger?

According to our simulations, the primary star, which is now an sdO star, will end its life as a CO white dwarf. The secondary star, which is now an ~10M_⊙ Be star, is sufficiently above the threshold where a single star ends its life in a core-collapse supernova. Without interaction, neither of the components would have been massive enough to produce a core-collapse supernova. As such, φ Persei is an example of a system that is expected to produce a late core-collapse supernova (Zapartas et al. 2017).

However, with a separation on the order of 200R_⊙, we expect the Be star (initially the secondary star) to fill its Roche lobe as it completes its main sequence evolution and swells up. This will initiate a phase of reverse Case B mass transfer from the secondary to the remnant of the primary, which is now a hot, massive young white dwarf. Since the donor star is considerably more massive than its companion, the orbit is expected to shrink rapidly in response to mass transfer. The resulting phase of mass transfer may thus be unstable, resulting in a common envelope phase. The outcome of this phase is uncertain. We expect that the most likely outcome is that the donor star engulfs the hot young white dwarf. The result is a merger with a peculiar interior composition. In Zapartas et al. (2017), we provide an extensive discussion and speculation of the possible outcomes of such a system.

5 Discussion and model variations

In this section, we discuss model variations, including variations in the mass transfer efficiency (Sect. 5.1) and different cases of mass transfer (Sect. 5.2), followed by alternative explanations for the excessive brightness of the sdO star in φ Persei (Sect. 5.3).

5.1 A case for conservative mass transfer

Our default models assume a mass transfer efficiency that is governed by the thermal rate of the accreting star, as explained in Sect. 3.1. The mass transfer efficiency is however one of the primary uncertainties in binary evolution (e.g., de Mink et al. 2007, and references therein). We therefore explore the impact of uncertainties in the mass transfer efficiency on the fit solutions. We explore model variations where we adopt different constant values for the mass transfer efficiency β ∈ {0.25, 0.5, 0.75, 1}. This affects the best-fitting solutions and confidence intervals that we find for the initial masses, initial orbital periods, and system age, as summarized in Table 2. The conservative models generally lead to solutions that have lower initial masses for both stars, larger initial orbital periods, and larger system ages.

Our default simulation and the simulations adopting a fixed value for β = 1 (conservative mass transfer) return identical solutions. Models adopting conservative mass transfer yield better-fitting solutions: we find ${\chi }_{\text{min}}^2=3.77$ for β = 1 (fully conservative mass transfer), while the fit solution is significantly worse for less conservative models with β = 0.25: ${\chi }_{\text{min}}^2=13.2$.

The preference for conservative mass transfer to explain the system parameters of φ Persei was pointed out earlier by Pols (2007). He used the mass of the sdO star to directly estimate the initial mass of its progenitor star assuming a moderate amount of core overshooting. The assumption that the progenitor of the sdO star was the initially more massive star of the system directly provides an upper limit to the initial mass of the secondary star. By comparing the inferred initial masses with the present day masses, Pols (2007) concluded that a mass transfer efficiency β > 0.7 is required, in agreement with our findings.

An independent argument in favor of conservative mass transfer comes from the inferred age. The φ Persei system is a likely member of the α Persei cluster, which has an estimated age of 52–60 Myr (see Sect. 4.2). This is consistent with the ages we obtain when assuming relatively conservative mass transfer. The best fit ages obtained with our non-conservative models, β ≤ 0.5, differ by more than 1σ from the age of the α Persei cluster. This further argues in favor of a relatively conservative mass transfer, but it will be interesting to see whether Gaia data can give a more definitive answer concerning cluster membership.

Finally, it is worth noting that our best fit models indicate a initial mass ratio of M_2,ini∕M_1,ini ≈ 0.5. Evolutionary simulations of interacting binaries typically report that the accreting star in such systems swells significantly during the mass transfer phase (Neo et al. 1977) such that it may possibly even fill its Roche-lobe as well. If true, this would imply that the system briefly evolved through a near-contact or contact configuration. It has been suggested that such systems would experience significant mass loss from the system or even result in a stellar merger (e.g., Pols 1994; Wellstein et al. 2001; de Mink et al. 2007). This is clearly in contradiction with the present-day appearance of the system. The current parameters of φ Persei provide a strong case that binary systems with similar initial mass ratios can evolve through a Roche lobe overflow phase without any significant loss of mass from the system.

5.2 Different evolutionary scenarios

The only solutions that we have found in our analysis to explain the present-day parameters of φ Persei all involve stable Case B mass transfer, i.e., mass transfer where the Roche lobe filling star is crossing the Hertzsprung gap (Kippenhahn & Weigert 1967).

Even when we extend our grid of initial periods to include much wider systems, we do not find additional solutions. Such wide systems undergo late Case B or Case C mass transfer, where the donor star is a giant with a convective envelope at the onset of mass transfer. Mass transfer is unstable for these systems and leads to a common envelope phase. The ejection of the envelope would lead to the orbit shrinking, which could in principle, with some fine-tuning, explain the present day orbital period. However, we find that such solutions are incompatible with the present day masses. The initial mass of the donor star is effectively fixed by the present-day mass of the sdO star, since it sets the core mass of the original primary star. This implies that the donor star initially cannot have been much more massive than about 7M_⊙. At the same time, solutions involving common envelope evolution imply that the secondary does not accrete a significant amount of mass. This means that the present-day mass of the secondary, which is now the Be star of 9.6 ± 0.3M_⊙, sets a lower limit to the initial mass of the primary. These two constraints are incompatible with each other.

We also do not find solutions resulting from shorter period systems that undergo Case A mass transfer (mass transfer from a main sequence donor). Such solutions imply an initial period of ~3.0 days or less, which are inconsistent with the large present-day orbital period of φ Persei (126 days) for any reasonable assumption of mass and angular momentum transfer. This can be shown directly using the analytical initial-final period relation (Eq. (16.20) in Tauris & van den Heuvel 2006, cf. Soberman et al. 1997), which requires as input the initial and final mass ratio and the mass transfer efficiency β. For angular momentum loss from the system, we again assume that all escaping matter takes away the specific angular momentum of the secondary star (see Sect. 3.1.1). We explore the parameter space with primary mass M_{1, 0} ∕ M_⊙∈{7.2, 13.0}, mass ratio q₀ ∈{0.0, 1.0}, and mass transfer efficiency β ∈ {0.0, 1.0}. For the current values of M₁ we use the stripped core mass which depends on M_{1, 0}. This stripped core mass is obtained by binary_c. The current M₂ follows from M_{1, 0}, q₀, and β. We find that the lowest possible initial orbital period for which P_{orb, obs} = 126 days and M_{2, obs} = 9.6M_⊙ can still be produced is P_{orb, 0} ≈ 4 days, which is larger than the maximum period for case A mass transfer at this mass. From these results, we conclude that we were not able to find solutions for a case A mass transfer scenario.

We have not considered possibilities involving a triple companion, but we see no reason why this would provide an explanation of the overluminosity of the stripped star.

5.3 Alternative explanations for the excessive luminosity of the stripped star

An alternative hypothesis that can be put forward to explain the high luminosity of the stripped star in the φ Persei SES (with respect to its mass) is that the stripped star is currently accreting matter. This is a possibility discussed by Mourard et al. (2015). Even though the visible circumstellar disk of the Be star does not extend to the vicinity of the SES in the case of φ Persei, this does not exclude the possibility of mass accretion onto the SES. We use a simple estimate to assess the likeliness of this hypothesis.

The stripped star appears to be almost an order of magnitude too luminous for helium core burning. Disk accretion would thus be the dominant energy source of the SES. Adopting its currently observed mass and radius, we estimate the accretion rate that is necessary to reach its observed luminosity using the following expression, $$L_{\text{acc},\text{\hspace{0.05em}SES}}=\frac{G{M}_{\text{SES}}{\dot{M}}_{\text{SES}}}{R_{\text{SES}}},$$(7)

and obtain an accretion rate of Ṁ > 10⁻⁴M_⊙ yr⁻¹. This is several orders of magnitude higher than the typical mass loss rate of Be stars (e.g., Waters 1986). We therefore deem this scenario implausible.

Another possible hypothesis is that the stripped star exhibits a remaining layer of hydrogen and that hydrogen shell burning could be responsible for the high inferred SES luminosity. A closer inspection of the detailed MESA models shows that the fractional hydrogen luminosity is around 1% at the onset of helium shell burning, and it decreases further afterwards. The total hydrogen mass in the object at the onset of helium shell burning was 0.004M_⊙. A much higherhydrogen mass than that would be needed for sustainable high luminosity as a result of hydrogen fusion. Götberg et al. (2017) show that even in the extreme case of very metal-poor stripped stars, which retain a significant layer of hydrogen, the relative contribution of hydrogen burning luminosity is at best comparable to that produce by helium burning. We consider it unlikely that this can explain an order of magnitude discrepancy in the luminosity. A further reason to deem this scenario unlikely is that a larger hydrogen burning envelope would bloat the star and reduce its effective temperature. The MESA model is already on the cold side of the observationally derived temperature.

6 Implications of the short-lived evolutionary stage of φ Persei for the population of undetected SESs

Given the large fraction of young stars with nearby binary companions that are expected to interact during their lives, we might naively expect SES systems to be relatively common. As a ball park estimate, we can assume that a fraction f of all stars of intermediate mass are members of a binary system that will produce an SES. We further safely assume that the duration of the SES phase is about 10% of the total lifetime, i.e., similar to the helium core burning lifetime. For any reasonable assumptions for f and the assumption of continuous star formation, this leads to the expectation that on the order of few percent of all stars of intermediate mass shouldharbor an SES (see Götberg et al. 2018). This ball park estimate appears to be in contrast with the very small number of systems that are currently known. In the remainder of this section we discuss the clues that the φ Persei evolutionary stage is giving us so that we can understand this discrepancy.

Fig. 5 Left panel: fraction of time spent by SESs in systems similar to φ Persei (see Sect. 4.4) in the helium core burning phase (pale blue) and the helium shell burning phase (bright blue). Right panel: fraction of time spent by SESs in systems similar to φ Persei below (pale green) and above (bright green) the luminosity limit Lobs. Before adding up all the time steps in the grid to obtain the plots above, each time step is weighted by the likelihood of fit from the derived initial parameters.

6.1 How rare is φ Persei?

In the previous sections we argue that the SES in the φ Persei system most likely resides in the short-lived evolutionary phase of helium shell burning. The brightness of the SES in φ Persei may have contributed to its detection and characterization as an SES in such a short-lived phase. For the full population of SESs, we expect the majority to reside in the long-lived helium core burning phase instead. Using our simulations we can estimate how many undetected systems we expect for each system like φ Persei.

We quantify this using our model grid by providing an average for the ensemble of systems that have birth parameters similar to those of φ Persei. To obtain this we compute the relative duration of helium shell burning for all SESs in our grid and weigh them with their likelihood of the fit to the derived initial parameters of φ Persei. We find that the helium shell burning phase lasts on average around 10% of the total helium burning lifetime (Fig. 5, left). During this phase the SES increases in brightness. In a similar fashion, we can compute the relative duration that SESs in systems with birth parameters similar to φ Persei are as bright as or brighter than the φ Persei SES. We find that an SES born in a system with similar birth parameters to φ Persei spends only ${2.4}_{-0.7}^{+0.9}\%$ of the SES lifetime having a luminosity L > L_obs (Fig. 5, right). The error bars here reflect the uncertainty in the observed luminosity. We obtain a very similar result if we analyze the MESA evolutionary track instead, which has a luminosity of L >L_obs for ${2.7}_{-1.3}^{+2.1}\%$ of its lifetime.

The stripped star in φ Persei thus belongsto the top ~2–3% brightest inits class. This implies that for each system similar to φ Persei, with a stripped star that is as bright or brighter than the stripped star in φ Persei, we expect that about 30–50 systems will host a dimmer stripped star, most of which reside in the helium core burning phase. Their presence will be very hard to detect through their luminosity alone. This, by itself, is providing insight in the observational difficulties to detect stripped stars.

6.2 The φ Persei system in the context of the known SESs

To try and place φ Persei in context we have conducted a literature study and have attempted to compile a complete overview of all currently known SESs with main sequence companions. We have only been able to find 11 systems matching these criteria¹. This includes φ Persei as well as the systems FY CMa and 59 Cyg, already used for comparison in Fig. 4. Unfortunately, in most cases several of the parameters are still unknown or only contradictory estimates are available. We provide the data that we were able to compile in Table 3.

We found a total of five systems that appear to be systems hosting a subdwarfs plus a Be star (Be+SdO). In addition to the three systems mentioned above, the recently detected SESs in 60 Cyg (Peters et al. 2016) and HD2142 (Wang et al. 2017) also belong to this group. They all host subdwarfs with masses around 1M_⊙, temperatures of 40–50 kK, high effective gravities of logg = 4.2–5. Their companions have masses ranging from 8–12M_⊙ and orbital periods ranging from 1 to 5 months. They all appear to be post-interaction systems resulting from stable and rather conservative mass transfer.

We found an additional set of five systems that have been marked to include helium giants (see Dudley & Jeffery (1993) and references therein, and Jeffery & Aznar Cuadrado 2001). These helium giants are thought to be helium stars that are in a very late phase of helium shell burning. The expectation is that these binaries are currently going through or are about to go through their second mass transfer phase (e.g., Schoenberner & Drilling 1983), as is predicted by De Greve & De Loore (1977) and Delgado & Thomas (1981). In this case, the donor star fills its Roche lobe again during a late helium shell burning phase in which it cools and expands. The SESs in these systems have lower temperatures, 9–30 kK, and lower effective gravities of log g = 1.5–3.5. We note that several parameters for these systems result from studies conducted in the 1980s. Often the nature of the companion is not known and/or its presence is only inferred from radial velocity variations. Observing these systems again would be worthwhile in order to obtain more reliable measurements. One of these stands out, υ Sgr, which is the most massive one in this category. The SES is a possible progenitor of a core-collapse supernova.

We have added one additional system to this compilation that does not belong to either of these groups. This is the quasi-WR star in HD 45166, which has a late B-type companion in a tight orbit of 1.6 days. Also worth mentioning is that Smith et al. (2018) and Götberg et al. (2017) recently argued that the new class of WN3/O3 stars discovered by Massey et al. (2014, 2015, 2017) may also be SES stars. At least one of these objects is a short-period eclipsing binary system (Graczyk et al. 2011). Whether all of these stars have companions is still a topic of debate.

It is striking that the majority of the SESs in these systems appear to be in their helium shell burning phase. This includes all five HeGB systems, φ Persei, and 59 Cyg. That is a total of at least seven out of a sample of eleven systems. For each of these seven we expect about 10–50 additional systems to still harbor undetected dimmer SESs. This seems to suggest we are missing hundred(s) of SESs in the sample from which these binary systems are drawn. The trend in recent observationsappears to agree with this interpretation. Only after a dedicated analysis did Peters et al. (2016) and Wang et al. (2017) find an sdO star with uncertain luminosity as a companion to a Be star. Moreover, Wang et al. (2018) found an additional 12 candidate Be+sdO systems, the stellar parameters of which remain to be determined.

At this stage the sample is far too heterogeneous and poorly characterized to derive any solid conclusions, and there may be further selection effects at play. It does however provide clear motivation for a more systematic search for these elusive SESs, especially given their astrophysical significance.

7 Summary and conclusions

A detailed study of nearby post-interaction binary systems is one of the important ways to provide crucial insight in the uncertain physical processes of binary interaction and to help us to better understand the population of the exotic products of binary interaction. Of particular interest is the apparent paradox of the scarcity of missing stripped stars. Even though binaries are common, very few systems are known that reside in the first long-lived post-mass transfer phase. The binary system φ Persei is thought to be one of these systems. It contains an ~1.2M_⊙ hot subdwarf and an ~9.6M_⊙ Be star in a 126.7-day orbit.

We systematically reanalyzed the binary system φ Persei using an extensive grid of binary evolutionary models and used this information to derive constraints on the progenitor, infer the current evolutionary state, and predict its future evolution. Our main findings are listed here:

• Our best fit for the progenitor of φ Persei is a post-interaction binary system where the hot subdwarf is the stripped core of a star with an initial mass of 7.2 ± 0.4M_⊙. For the rapidly rotating Be star we derive an initial mass of 3.8 ± 0.4M_⊙. We derive an initial orbital period of 16 ± 4 days. We find that models including conservative mass transfer, i.e., models where at least 75% of the mass lost by the donor is accreted by the secondary, are favored, consistent with earlier findings.

• We derive an age of 57 ± 9 Myr for the system. This is in excellent agreement with the most recently determined age of its suspected host cluster α Persei, 52–60 Myr. This would imply that the 9.6 ± 0.3M_⊙ Be star in φ Persei is a rejuvenated blue straggler, which has a turn-off mass of about 7M_⊙. This also makes φ Persei a likely progenitor of a late core-collapse supernova, i.e., one that explodes after all massive single stars in the cluster have already ended their lives.

• We find that the subdwarf is almost an order of magnitude brighter than expected for a helium core burning star of this mass. It is also ~3 times more luminous than the subdwarf with a similar mass in the system FY CMa. We can only reproduce the excessive luminosity with our models if the subdwarf has already completed helium core burning and is now in a more luminous phase where it is burning helium in a shell around the core.

• The surprisingly large luminosity of the subdwarf in φ Persei seems to imply that we have caught the system in a short-lived late evolutionary stage. Based upon lifetime arguments, we expect only 2–3% of subdwarfs of this mass to be this bright. This would imply that there are up to 50 systems with birth parameters similar to those of φ Persei that contain a less evolved subdwarf which has remained undetected due to various observational biases. This appears to be an important clue for understanding the apparent paradox of the scarcity of detections of stripped stars in binary systems.

• We compiled an overview of all currently reported subdwarfs and helium giants with main sequence companions in literature. The sample only yields 11 systems. At least seven of these systems appear to host stripped stars that currently reside in the late stages of the short-lived helium shell burning phase. If true, this would also indicate that we are still overlooking the presence of possibly hundred(s) of less evolved stripped stars.

Understanding the apparent paradox of the scarcity of detections of stripped stars in binary systems is important for various fields in astrophysics. It has been suggested that stripped stars play an important role by contributingfar-UV and ionizing radiation in older stellar populations. Higher mass counterparts of the φ Persei stripped star are considered to be Ib∕c supernova progenitors. As such, binary systems that eventually produce gravitational wave sources through the classic formation scenario are expected to evolve through a post-interaction phase very similar to that occupied by φ Persei.

Our results imply that there is a substantial population of systems analogous to φ Persei that containless evolved stripped stars and that have evaded detection so far. Dedicated searches for these systems will be extremely valuable since they can provide strong constraints on the physics of binary interaction, as we have shown.

Acknowledgements

The authors would like to thank Rob Izzard for providing the binary_c code; Ulrich Heber, Pablo Marchant, and Philipp Podsiadlowski for help and input; and Norbert Langer for useful comments on the draft. SdM acknowledges support from a Marie Sklodowska-Curie Action Incoming Fellowship (H2020 MSCA-IF-2014, project BinCosmos, id 661502). The authors further acknowledge the Leiden Lorentz Center for support of the workshop “The Impact of Massive Binaries Throughout the Universe”.

References

All Tables

All Figures

Fig. 1 Posterior probability distributions of the fit parameters (yellow histograms) of the φ Persei system together with a fitted Gaussian curve (in black). For comparison, the observed parameters for the φ Persei system (bluevertical dashed lines) and their 1 σ uncertaintyintervals (blue shaded regions) are shown. Top panels: current masses of the sdO and Be star, MSES and MBe. Central panels: effective temperatures, Teff, SES and Teff, Be. Bottom panels: present-day orbital period Porb and the ratio of their radii RSES∕RBe.

In the text

	Fig. 2 Posterior probability distributions of the birth parameters (green histograms) of the φ Persei system together with a fitted Gaussian curve (in black). Top panels: birth masses of the sdO star and the Be star, MSES, 0 and MBe, 0. Bottom panels: initial orbital period Porb, 0 and the system age.
In the text

	Fig. 3 Hertzsprung–Russell diagram showing the evolutionary track of a stripped envelope star from the best-fitting system derived with binary_c. This object has a mass of 1.36M⊙. It is simulated with both binary_c (dashed line) and MESA (solid line). Models in the helium core burning phase are displayed in blue, models in the helium shell burning phase in orange.
In the text

Fig. 4 Evolution of the luminosity and mass of stripped envelope stars of 0.8M⊙, 1.2M⊙, and 1.6M⊙. Helium core burning models are displayed in blue, while helium shell burning models are displayed in orange. The blue shaded area indicates the parameter space in which the models predict SESs to be helium core burning; in the area above, SESs are predicted to be helium shell burning. In the gray hatched region, the models do not predict helium burning SESs. Top (solid lines): MESA tracks. Bottom (dashed lines): binary_c tracks.

In the text

Fig. 5 Left panel: fraction of time spent by SESs in systems similar to φ Persei (see Sect. 4.4) in the helium core burning phase (pale blue) and the helium shell burning phase (bright blue). Right panel: fraction of time spent by SESs in systems similar to φ Persei below (pale green) and above (bright green) the luminosity limit Lobs. Before adding up all the time steps in the grid to obtain the plots above, each time step is weighted by the likelihood of fit from the derived initial parameters.

In the text