A Be-type star with a black hole companion

Stellar-mass black holes have all been discovered through X-ray emission, which arises from the accretion of gas from their binary companions (this gas is either stripped from low-mass stars or supplied as winds from massive ones). Binary evolution models also predict the existence of black holes accreting from the equatorial envelope of rapidly spinning Be-type stars (stars of the Be type are hot blue irregular variables showing characteristic spectral emission lines of hydrogen). Of the ~80 Be X-ray binaries known in the Galaxy, however, only pulsating neutron stars have been found as companions. A black hole was formally allowed as a solution for the companion to the Be star MWC 656 (also known as HD 215227), although that was based on a single radial velocity curve of the Be star, a mistaken spectral classification and rough estimates of the inclination angle. Here we report observations of an accretion disk line mirroring the orbit of the Be star. This, together with an improved radial velocity curve of the Be star through fitting sharp Fe II profiles from the equatorial disk, and a refined Be classification (to that of a B1.5-B2 III star), reveals a black hole of 3.8 to 6.9 solar masses orbiting MWC 656, the candidate counterpart of the gamma-ray source AGL J2241+4454. The black hole is X-ray quiescent and fed by a radiatively inefficient accretion flow giving a luminosity less than 1.6 x 10-7 times the Eddington luminosity. This implies that Be binaries with black-hole companions are difficult to detect by conventional X-ray surveys.

displays considerable scatter, likely to be caused by filled-in emission from the circumstellar wind contaminating the broad absorption profiles (a common limitation in the analysis of BeXBs: see, for example, ref. 10) and only a tentative orbital solution is available 5 . Here we revisit the 32 Liverpool telescope spectra previously reported 5 . These are complemented with 4 additional Liverpool telescope spectra and a further high-resolution echelle spectrum obtained with the 1.2-m Mercator telescope (see Methods and Extended Data Table 1).
A close-up of the Mercator telescope spectrum is presented in Fig. 1, showing classic Fe II emission lines from the Be circumstellar disk. In addition, a He II 4,686 Å emission line (which was overlooked in a previous work 5 ) stands out clearly: its presence is remarkable because it requires temperatures hotter than can be achieved in disks around B-type stars. Further, the He II profile is double-peaked, which is the signature of gas orbiting in a Keplerian geometry 11 .
Gaussian fits to the He II profiles in the Liverpool telescope spectra reveal that the centroid of the line is modulated with the 60.37-day orbital period, reaching maximum velocity at photometric phase 0.06 (see Methods and Extended Data Fig. 1). This is approximately in antiphase with the radial velocity curve of the Be star 5 , a strong indication that the He II emission arises from gas in an accretion disk around the invisible companion and not from the Be disk. We can therefore use its radial velocity curve to trace the orbit of the Be companion.
An eccentric orbital fit to the He II velocities was performed using the Spectroscopic Binary Orbit Program (SBOP 12 ), fixing the period to 60.37 days (Methods); the resulting orbital elements are given in Extended Data Table 2. The orbital evolution of the He II line is presented in Fig. 2. The line flux is also found to be modulated with the orbital period (Methods and Extended Data Fig. 1), owing to the presence of an S-wave component swinging between the double peak (see Fig. 2).
In order to improve on the radial velocity curve of the Be star previously reported 5 , we fitted the sharp double-peaked profile of the Fe II 4,583 Å emission line with a two-Gaussian model (Methods). Fe II lines are known to arise from the innermost regions of the circumstellar disk 13,14 , and therefore reflect the motion of the Be star much more accurately than the broad He I absorptions. The Fe II velocities were also modelled with SBOP, fixing the period to 60.37 days. The eccentric orbital fit results in orbital elements that are consistent with the Fe II orbit being the reflex of the He II orbit, as expected from the motion of two components in a binary system (see Methods and Extended Data Table 2). Only the eccentricities are slightly discrepant, but just at the 1.6σ level.
Consequently, we modelled the ensemble of Fe II and He II radial velocities with a doubleline eccentric binary orbit in SBOP. Figure 3 presents the radial velocity curves of the two emission lines with the best combined solution superimposed. The resulting orbital elements are listed in Table 1. Our solution yields a mass ratio q = M 2 /M 1 = 0.41 ± 0.07, which implies a rather massive companion star. A precise determination of the companion's mass requires an accurate spectral classification of the Be star. Accordingly, we have compared the Mercator telescope spectrum with a collection of observed B-type templates, broadened by 330 km s −1 to mimic the large rotation velocity in MWC 656 (Methods and Extended Data Fig. 3). Using excitation temperature diagnostics based on several absorption line ratios, we determine a spectral type of B1.5-B2, and the strength of the metallic lines combined with the moderate width of the Balmer absorption wings implies luminosity class III (Methods). Our adopted B1.5-B2 III classification implies a mass of 10-16 solar masses (M Í ) for the Be star (Methods), and hence a companion star of 3.8-6.9 M Í .
The large dynamical mass of the companion to the Be star MWC 656 is puzzling. A normal main sequence star with such mass would have a spectral type in the range B3-B9 and its spectrum would be easily detected in the optical range. Nor can it be a subdwarf, because these typically have masses in the range 0.8-1.3 M Í (ref. 15). The stripped He core of a massive progenitor (that is. a Wolf-Rayet star) is also rejected because it possesses strong winds which show up through intense high excitation emission lines, not present in our spectra. In addition, this should be detected as an ultraviolet excess in the spectral energy distribution, which is not observed in the fluxes available in the literature (Methods). On the other hand, the evidence of a He II accretion disk encircling the companion star strongly points towards the presence of a compact object. The large dynamical mass rules out a white dwarf or a neutron star, so the only viable alternative is a black hole. It should be noted that none of the ~170 BeXBs curently known 4 shows any evidence for an accretion disk, providing circumstantial evidence for a difference in the nature of the compact stars. The accretion disk in MWC 656 is expected to also radiate Balmer and He I lines but these are blended with the corresponding (stronger) emission lines from the Be disk and thus are not detected.
MWC 656 is a key system in the study of BeXBs and massive binary evolution. At a distance d = 2.6 ± 0.6 kpc (Methods) it is relatively nearby and also one of the visually brightest Be binaries 7 . It seems thus reasonable to assume that many other Be-BHs exist in the Galaxy but remain hidden by the lack of transient X-ray activity. Analysis of archival ROSAT images yields an upper limit to the X-ray flux at energies of 0.1-2.4 keV of 1.2 × 10 −13 erg cm −2 s −1 (Methods) which, for our estimated distance, translates into an X-ray luminosity L x < 1.0 × 10 32 erg s −1 or < 1.6 × 10 −7 times de Eddington luminosity, L Edd . Therefore, accretion is highly inefficient in MWC 656, akin to accretion on black holes in quiescent low-mass X-ray binaries 16 , where accretion disks are truncated at ~10 2 -10 4 Schwarzschild radii and then behave as an advection dominated accretion flow 17 .
In the context of disk instability theory, the very low mass-transfer rates expected for BeXBs (with peak values of ~10 −11 M Í yr −1 near periastron) lead to extremely long outburst recurrence periods or even to completely suppressed transient activity 18 . It is the dormant condition of the accretion disk together the absence of a solid surface reradiating the accretion energy that makes Be-BHs very difficult to detect through X-ray surveys, thus providing an explanation for the missing Be-BH population. This is in stark contrast with the other black-hole high-mass X-ray binary known in the Galaxy, Cygnus X-1, where an X-ray persistent accretion disk is fed by the powerful wind (~ 10 −8 M Í yr −1 ) of an O supergiant star 19 .
The detection of a Be-BH is also important to our understanding of BeXB evolution.
Whereas the total number of neutron-star BeXBs in the Galaxy depends strongly on the distribution of kick velocities, the number of Be-BHs is very sensitive to the survival probability during the common envelope phase 3 . Modern population synthesis models predict a Galactic number ratio of neutron-star to black-hole BeXBs of 54, for the case of no common envelope survival during the Hertzsprung gap and a Maxwellian distribution of kick velocities with reduced root mean square σ =133 km s −1 (model C in ref. 3). There are currently ~81 BeXBs known in the Galaxy with ~48 pulsating neutron stars 4,20 , and thus our discovery of a black-hole companion to MWC 656 is consistent with these model predictions. However, it should be noted that the X-ray spectra of the remaining BeXBs, whenever they are available, also indicate the presence of a neutron star. Further, in stark contrast with the known BeXBs, MWC 656 has been identified through a claimed γ-ray flare (see Methods) and not by its X-ray activity. This seems to imply that the discovery of Be-BHs is observationally biased, in which case common envelope mergers would be less frequent than commonly assumed and/or neutron star kicks would be best described by the radio pulsar birth velocity distribution 3 . Last, it is interesting to note that MWC 656 will probably evolve into a black-hole/neutron-star binary, a potential source of strong gravitational waves and a short γ-ray burst (Methods).

METHODS SUMMARY
The Fe II line was fitted with a two-Gaussian model with Gaussian positions, intensities and separation left as free parameters. The Fe II velocities were obtained from the mean of the Gaussians offset with respect to the rest velocity at 4,583.837 Å. A detailed radial velocity study of the He II profile was performed using the double-Gaussian technique 21 , and this shows that the systemic velocity is pushed down in the line core by the S-wave component while the phasing and velocity semiamplitude remain very stable. Therefore, a +30 km s −1 offset was applied to the He II velocities before fitting all the data points with a double-line orbital model.
The S-wave also modulates the He II flux with the orbital period, and its phasing can be interpreted as either enhanced mass transfer during periastron or the visibility of a hotspot in the outer accretion disk. The spectral classification of the star was obtained by direct comparison with a range of templates conveniently broadened, and the best match is provided by the B1.5 III star HD 214993. We used several calibrations available in the literature to constrain the mass of MWC 656 from its spectral type, including evolutionary tracks, a dynamical determinations from the detached eclipsing binary V380 Cyg and robust lower limits from dynamical masses of main sequence stars. The distance is obtained through combining the absolute magnitude of B1.5-2 III stars with the observed brightness of MWC 656, corrected for intestellar reddening.
Upper limits to the X-ray flux of MWC 656 are derived from archival ROSAT and Swift pointings, using a neutral hydrogen column density N H = 1.4 × 10 21 cm −2 and a photon index Γ= 2.0, typical of black holes in quiescence. We futher discuss the future evolution of MWC 656 and its fate, a possible black-hole/neutron-star binary. The solution was obtained from a combined fit to the radial velocity curves of the He II 4,686 Å and Fe II 4,583 Å lines. The orbital period P orb has been fixed to the photometric value 6 Table 1.

Radial velocity analysis
Radial velocities were extracted from the He II 4,686 Å emission by fitting a single Gaussian to the line profile in the 36 Liverpool telescope spectra. A least-squares sine fit to the velocity points yields an orbital period of 59.5 ± 0.6 days (Extended Data Fig. 1a) which is consistent within 1.4 σ with the (more accurate) photometric period determination 6 . The difference is also explained by the fact that our radial velocities only cover two full orbital cycles. Consequently, we henceforth adopt 60.37 ± 0.04 days as the true orbital period of the binary. The radial velocities were modelled with an eccentric binary orbit using the code SBOP 12 , with the period fixed to 60.37 days. Individual points were weighted proportionally to 1/σ 2 , where σ is the radial velocity uncertainty. We adjusted the following orbital parameters: epoch of periastron (T 0 ), eccentricity (e), argument of the periastron (ω), systemic velocity (γ) and velocity semiamplitude (K). The resulting orbital elements are listed in the first column of Extended Data Table 2, together with their implied fundamental binary parameters. The fitted solution yields maximum velocity at phase 0.06, where phase 0 is arbitrarily set to HJD 2453243.3 or the epoch of maximum brightness in the photometric light curve 6 . This is approximately in antiphase with the radial velocity of the Be star previously determined 5 , and therefore the He II emission most probably follows the orbit of the companion star.
To further constrain the orbital motion of the Be star we measured radial velocities from the Fe II 4,583 Å emission line by fitting a two-Gaussian model to the individual profiles. The Gaussians are set to have identical width but their intensities and separation are allowed to vary in order to account for possible profile changes. In any case, these three free parameters are found to be very stable, with differences which are always within 10 %. Only in one case was the double-peak separation found to be 13% smaller than the average. This corresponds to one spectrum where the profile appears blurred. For this particular case we decided to fix the Gaussian separation to the average value, 270 km s −1 . The radial velocity of the Fe II line was then obtained from the center of the double Gaussian model relative to the line rest velocity at 4,583.837 Å. We subsequently fitted the radial velocity points with an eccentric orbit model using SBOP after fixing the orbital period to 60.37 days and weighting data points proportionally to 1/σ 2 , as before. The best set of fitting parameters are listed in the second column of Extende Data Table 2.
The two independent solutions presented in Extended Data Table 2 give periastron phases These are all strong indications that both radial velocity curves display the reflex motion of two components in a binary system, thus endorsing fitting a combined double orbit to the ensemble of 72 He II and Fe II data points. In any case, the radial velocity amplitudes of the single line orbits are found to be very robust and in excellent agreement with the combined double-lined solution presented in the main text. These are the key parameters constraining the binary mass ratio and thus the nature of the companion to the Be star.

The diagnostic diagram
Extended Data Table 2 shows a 30 km s −1 difference between the systemic velocities of the Fe II and He II solutions, which we interpret as due to contamination by the S-wave component seen in the core of the He II profile 22 . To test this, we measured radial velocities from different velocity bands of the He II line profile, using the double-Gaussian technique 21 . To further increase the signal-to-noise ratio, we averaged the 36 Liverpool telescope spectra into 15 phase bins, relative to the photometric ephemeris. Every phase-binned spectrum was convolved with a two-Gaussian passband where the Gaussian separation was varied between a =200-600 km s −1 in steps of 50 km s −1 . Gaussian separations a < 200 km s −1 produce radial velocities which are corrupted by the variable shape of double-peaked profile, while for a > 600 km s −1 the line flux becomes too weak. Each radial velocity curve was subsequently fitted with the expression V(φ) =γ+K sin 2π(φ−φ 0 ), where φ is the binary phase, and the evolution of the fitting parameters K, γ and φ 0 versus the Gaussian separation a is displayed in Extended Data Fig. 2 as a diagnostic diagram 23 . Note that we prefer to fit a simple sine wave model rather than a full eccentric orbit because the extra fitting parameters (that is, e, ω and periastron phase) become poorly constrained as we approach the noisy wings of the line profile. Given the small orbital eccentricity the fitted sine wave parameters are still meaningful, although they should be taken as a first approximation to the true values because of the model oversimplification.
The high velocity wings of the He II profile are formed in the inner regions of the accretion disk and are less affected by the core S-wave component. Therefore, they are expected to closely trace the motion of the compact star. Extended Data Fig. 2 shows that, as we move away from the line core, the systemic velocity rises quickly from about −46 km s −1 to about −25 km s −1 , thus approaching the value of the Fe II solution. This demonstrates that the Fe II velocities provide a more accurate description of the true binary systemic velocity than the centroid of the He II line and, therefore, we decide to apply a +30 km s −1 offset to the latter. Beyond a > 500 km s −1 the continuum noise begins to dominate the radial velocities, as indicated by the steeper rise in the control parameter σ(K)/K (ref. 23), and thus the fitted parameters are less reliable. The diagram also illustrates that both phasing and velocity semiamplitude are very stable in the interval a =200-500 km s −1 , ranging between φ 0 =0.77-0.82 and K =76-92 km s −1 . Beyond a > 500 km s −1 K drops below 70 km s −1 and hence there is a possibility that the velocity semiamplitude of the compact star was overestimated in Table 1. Note, however, that this would raise the binary mass ratio which would make even stronger the argument for a black hole companion.

Orbital modulation of the He II Flux
Extended Data Fig. 1b shows that the equivalent width of the He II 4,686 Å line is strongly modulated with the orbital period. It peaks at phase ~0.9, which is very close to the maximum in the photometric light curve (that is, phase 0 by convention). On the other hand, the amplitude of the equivalent width modulation is an order of magnitude larger (~40 %) than that of the photometric light curve 6,24 . These two arguments imply that the equivalent width variability is driven by true changes in the line flux rather than in the continuum.
Extended Data Fig. 1 also reveals that the maximum equivalent width (phase 0.93 ± 0.04) almost coincides with the peak in the He II radial velocity curve, when the compact star is receding from us at maximum speed. The latter agrees well with the time of maximum visibility of the hotspot or shock region between the gas stream and the accretion disk. Alternatively, the modulation of the He II flux can be interpreted as caused by enhanced mass transfer near periastron, which is constrained to phase 0.01 ± 0.10 by our orbital solution.

Spectral classification, mass and distance to MWC 656
The spectral classification of MWC 656 is complicated by the effects of fast rotation, and by the presence of many (mostly Fe II) emission lines affecting many of the features useful for this purpose. To provide an improved classification, our high-quality Mercator telescope+HERMES spectrum was compared to the spectra of several MK standards taken with the same instrumentation and set-up as for MWC 656. Before this, we measured the rotational broadening in MWC 656 by applying the Fourier technique 25 , combined with a goodness-of-fit method 26 Fig. 3, the best match to the overall spectrum is provided by the B1.5 III standard HD 214993, although MWC 656 has rather stronger N II features. Nitrogen enhancement is frequently found in fast rotators 28,29 , and understood as a natural product of stellar evolution 30 . Nitrogen enhancement is often accompanied by C depletion. Therefore, we cannot completely rule out the possibility that MWC 656 is a B2 III giant with some C depletion (compare to the spectrum of the MK standard HD 35468 in Extended Data Fig. 3), but the strength of Si III and O II lines, which are the features most sensitive to temperature in this spectral range, strongly supports a B1.5 III classification.
Armed with the spectral classification we can now set constraints to the mass of the Be star However, the modulation has a semiamplitude of only 0.020 mag, suggesting that these other contributions have a very modest effect 6,24 and, therefore, we can safely conclude that the Be star dominates the observed optical flux. In any case, the large error in our distance estimate encompasses the uncertainty introduced by the small contribution of non-stellar light sources to the observed optical flux. Incidentally, the SED does not show any ultraviolet excess as would be expected if the companion to the Be star were a stripped He core. This gives additional support for the presence of a black hole in MWC 656.

Gamma-ray association
MWC 656 has been proposed as the optical counterpart of the transient GeV γ-ray source AGL J2241+4454 (refs 5,6). The association is, however, uncertain because AGL J2241+4454 was only detected during a 2-day activity period and it has a position error circle radius of 0.6 o (ref. 9). A few binary systems have been detected at GeV and/or TeV energies 43 , all showing orbitally modulated γ-ray emission and most of them thought to contain young non-accreting pulsars. In contrast, the accreting black hole Cygnus X-1 showed evidence of TeV emission during less than one day 44 , and also two different short transient episodes of GeV emission 45 .
The black-hole nature of MWC 656 and its putative flaring γ-ray activity are reminiscent to those of the well known black hole Cygnus X-1 but, remarkably, the accretion luminosity in MWC 656 is typically lower by more than ~5 orders of magnitude.
A candidate black-hole/neutron-star progenitor The future evolution of MWC 656 will probably lead to a black-hole/neutron-star binary 46 .
During the red giant phase the 13 M Í Be star will expand by several hundred solar radii 47 , thus engulfing the black hole. Mass transfer from the expanding Be star onto the black hole will be dynamically unstable and a common envelope will ensue. This is a highly dissipative process which leads to spiral-in of the black hole, efficient circularization of the orbit and the ejection of the Be star envelope. The outcome of the common envelope phase will then be a 2.9 M Í He star 48 and the present ~5 M Í black hole companion in a close circular orbit. In the event of a symmetric core collapse the new born neutron-star/black-hole binary will remain bound because less than half the total initial mass is expelled in the explosion 49 . In the case of an asymmetric supernova, the binary survival will depend on the magnitude and direction of the kick.
Black-hole/neutron-star binaries, which have not been detected yet, are instrumental in providing fundamental tests of gravitational theories, strong sources of gravitational waves and prime candidates for the production of short γ-ray burst through coalescence [50][51][52] . The fate of MWC 656 as a possible black-hole/neutron-star binary is very relevant because it provides tight empirical constraints on detection rates for gravitational wave observatories, such as advanced LIGO/VIRGO 53 .