Supersoft X-Ray Sources Identified with Be Binaries in the Magellanic Clouds^*

Very luminous supersoft X-ray sources (SSSs) have been observed since the end of the 1970s with the Einstein satellite, but they still pose unsolved riddles. More than half of the SSSs are transient sources, and a majority turn out to be post-outburst novae, in which the white dwarf (WD) keeps on burning hydrogen for a period of time ranging from a week to years after the outburst, with an atmospheric temperature of up to a million K. Novae are observed as SSSs in the Galaxy and the Magellanic Clouds, and now routinely even in M31 (for lists and reviews of novae in the Galaxy and the Local Group, see, among others, Orio et al. 2010; Henze et al. 2011, 2014; Orio 2012; Osborne 2015). However, the nature of numerous SSSs is not yet understood and may hold the key to outstanding astrophysical problems. Because of the large intrinsic luminosity, these sources are observed in the direction of external galaxies, in the Local Group and beyond, up to a distance of 15 Mpc (e.g., Di Stefano & Kong 2003, 2004) in regions of the sky affected by low absorption. Many non-nova SSSs are likely to be in close binaries hosting the hottest, most massive accreting and hydrogen-burning WDs, possibly on the verge of Type Ia supernova (SN Ia) explosions (see, for instance, van den Heuvel et al. 1992). The phenomenology of SSSs has an overlap with accreting black holes (see, e.g., Liu & Di Stefano 2008), and the most luminous SSSs may host, instead of WDs, massive stellar black holes accreting at super-Eddington luminosity.

A group of SSSs appears particularly interesting, potentially explaining the dependence of the SN Ia star formation rate (e.g., Sullivan et al. 2006) even without invoking double degenerate systems: the SSSs that are associated with massive stars, proposed to be Be-type stars (see review by Orio 2013). In addition to four SSSs whose new observation we describe in this paper, there is a fifth SSS in the Magellanic Clouds: the transient source RX J0527.8–6954 (Greiner et al. 1991, 1996), whose X-ray light curve is like that of a classical nova, has been positively identified with a Be optical counterpart (Oliveira et al. 2010). This source is in a very crowded field, and other possible counterparts were previously proposed, but the identification by Oliveira et al. (2010) with a B5e V star with subarcsecond and possibly bipolar Hα emission appears very convincing. This object also shows weak [O iii] emission. Moreover, among almost 100 SSSs in M31, in a statistically significant fraction of them have only massive young stars as detectable candidate optical counterparts (Orio et al. 2010).

Be stars with WD companions have not been detected yet, but there is reasonable circumstantial evidence that the rapidly rotating B star Regulus, for instance, hosts a WD companion (Gies et al. 2008). The mass range of B stars spans from a little above 2 to 20 M_⊙, so most of these stars end as WDs. Thus, the binary may have progenitor components in this range and does not require an unusual evolutionary path with large mass loss or mass transfer to end as a WD+Be system. However, mass transfer does also allow for unusual configurations in which the initial secondary, or less massive star, ends its life as a WD before the more massive companion. Matson et al. (2015) described observational evidence obtained for KOI-81, a system in which the original donor star was completely stripped and the secondary (a B star) was spun up. Wang et al. (2018) discovered many cases of Be stars with hot, stripped-down subdwarfs detected in the ultraviolet.

Waters et al. (1989) predicted a large fraction of Be+WD systems, and Raguzova (2001) calculated that 70% of Be systems have a WD companion, with most orbital periods less than a year. Such systems are not long-lived as SSSs but altogether last longer than neutron star+Be star systems, of which many are known (about 90 are known in the SMC and 33 in the LMC; Maravelias et al. 2014; Antoniou & Zezas 2016). McSwain & Gies (2005) found that more than 70% of Be binaries are spun up in the course of evolution and are expected to end as close binaries. Only one Be star with a black hole companion has been discovered (Casares et al. 2014), but it does not emit luminous supersoft X-rays; it is instead a rather hard X-ray source of low X-ray luminosity (1.6 × 10⁻⁷ times the Eddington luminosity).

We obtained low-resolution spectra of the proposed optical counterparts of four SSSs in the Magellanic Clouds with the Southern African Large Telescope (SALT; Buckley et al. 2006) 10 m telescope, and we show in this paper that all four have characteristics of Be stars. Two of the sources had already been observed spectroscopically, so we could assess whether the spectra evolved. In Section 2, we present the targets. The observations and reduction procedures are described in Section 3. In Section 4, we show the analysis, and we discuss the results in Section 5. A summary and conclusions are presented in Section 6.

None of our targets were persistent X-ray sources. Three of them were transients but have not been observed often enough to rule out recurrence. We know that one, XMMU J010147.5–715550, is a variable, recurrent SSS.

2.1. Previously Proposed Identifications without Optical Spectra

2.2. Proposed Identification with Published Optical Spectroscopy

We observed the four proposed optical counterparts described above between 2016 September and October with the Robert Stobie Spectrograph (RSS; Burgh et al. 2003; Kobulnicky et al. 2003) and between 2017 June and October with the echelle High Resolution Spectrograph (HRS; Bramall et al. 2010, 2012; Crause et al. 2014) mounted at the SALT.

The observations are summarized in Table 1. XMMU J010147.5–715550, SUZAKU J0105–72, and XMMU J052016–692505 were observed with the PG0900 grating in long-slit mode with a slit field of view (f.o.v.) of 1'' × 8^'. The resulting spectral resolution is about R = 1100 with a dispersion of about 0.97 Å pixel⁻¹, and the wavelength range is about 4060–7130 Å. MAXI J0158–744 and XMMU J052016–692505 were observed with the PG2300 grating in long-slit mode with the same slit. The resulting spectral resolution is about R = 2900 with a dispersion of 0.32 Å pixel⁻¹. The wavelength range is 4440–5460 Å. All of the targets were observed with the echelle HRS with its low-resolution grating, which allows us to observe the 3700–8870 Å spectral range with spectral resolution R = 15,000 and a dispersion between 0.024 and 0.045 Å pixel⁻¹.

The extraction, reduction, and analysis of spectra were performed with the Image Reduction and Analysis Facility (IRAF; version 2.16.1⁵ ). The spectra provided by SALT are already bias-subtracted. We extracted the RSS spectra with the apall task, wavelength-calibrated and flux-calibrated using the spectrophotometric standard stars LTT6248, EG274, LTT7379, and HILT600 and the usual IRAF reduction tasks. The standard stars were observed with the widest slit (width = 4'') and in the same position along the slit as the targets. We reduced the echelle spectra using the echelle package in IRAF. We performed the tracing using the flat-field frames, which have higher signal-to-noise ratios (S/Ns) than the objects. The objects were corrected for cosmic rays and extracted. Then, wavelength calibration was applied using a thorium–argon lamp, and the sky contribution for each aperture was subtracted. The spectra were flux-calibrated using the HR 7596, HR 5501, and HR 14943 standard stars and the usual standard and sensfunc tasks.

Absolute flux calibration cannot be achieved with the SALT telescope because it has a variable pupil and the illuminating beam changes during the observations; therefore, we only obtained relative flux calibration. The spectra shown in this work are either flux-normalized spectra using the flux values at λ = 5500 Å (RSS with PG0900) or λ = 4950 Å (RSS with PG2300) or flux-calibrated spectra (HRS), where the calibration is relative (not absolute).

In the RSS spectra taken with the PG2300 grating and in the HRS spectra, we observed double-peaked line profiles, which are typical of the rotating disks of the Be stars (see, among others, recent work by Labadie-Bartz et al. 2018). In the RSS spectra taken with the PG0900 grating, we could not resolve two peaks, even if they were present, due to low resolution. In all RSS spectra, we measured the velocity of the centroid of the lines by means of the IRAF task emsao. Because the line profiles are complex in the HRS spectra, we could not apply emsao. We estimated the radial velocity by measuring the position of both the blue/violet (V) and red (R) peaks; then, we took the average of the two peaks. We analyzed all of the spectra with the IRAF task splot to measure the equivalent width (EW) and fit the lines. We fitted the emission lines with one or more Gaussian functions and the underlying continuum with a straight line. To estimate the errors, we chose 10 different continuum levels to measure EWs, and we fitted the profiles 10 times, thus obtaining the mean and root mean square (rms) values for the EW, FWHM, and centroid positions. For double-peaked profiles, we also measured V/R, defined as the ratio of the V and R peak intensity (to compare our work with that of Dachs et al. 1992). The mean EW and rms values are given in Table 2, while kinematics data and V/R are in Table 3.

We fitted the lines that were clearly split, with an R and a V peak (PG2300 RSS and HRS spectra), with two Gaussian functions to measure the velocity difference (Δv) of the V and R peaks. Presumably, assuming the lines arise in a disk, these velocity differences represent 2v_d sin i, where v_d is the Keplerian velocity of a characteristic radius in the outer disk where the line flux is a maximum, and i is the disk inclination to the line of sight. When two peaks were not clearly measurable, for all PG0900 RSS spectra and some lines in the HRS spectra, we assumed that an estimate for Δv was given by the velocity dispersion value σ = 0.424 × FWHM. The October spectra of XMMU J052016–692505 obtained with the HRS show complex lines in which three peaks are detected, so we measured the FWHM directly on the line profiles and derived σ from the FWHM. In Figure 1, we show examples of how we fitted the lines observed with the different gratings or spectrographs. The wavelength-calibration errors were calculated by assuming a mean error of 20% of the dispersion value: the mean errors are 11, 4, and 2 km s⁻¹ for the PG0900, PG2300, and HRS data, respectively.

Zoom In Zoom Out Reset image size

4.1. Analysis of the Spectra

All of our RSS spectra show a steep blue continuum with Balmer hydrogen emission lines (Hα, Hβ, and Hγ), as expected for Be stars; in some cases, we also detected He i (λ = 5876, 6678, and 7065 Å) and Fe ii (λ = 4253, 5018, and 5317 Å) emission lines, and some hydrogen emission lines appear to arise from shallower, broader absorption features. In some of the spectra, the following hydrogen and helium absorption lines are also present: Hδ, He i (λ = 4143, 4168, 4387, 4471, and 4921 Å), and He ii (λ = 4200 Å). All of the RSS flux-normalized spectra are plotted in Figure 2 (PG0900 grating) and Figure 3 (PG2300 grating) to show which emission and absorption lines are present and measurable in which spectra. In both cases, an arbitrary constant was added to the flux values to make the comparison easier. The values of the normalization flux are indicated in the captions of the figures. We listed in Table 2 the measured EWs for the emission lines of each spectrum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Two spectra of XMMU J010147.5–715550 obtained 22 days apart clearly confirm the Be classification. In Figure 2, it is evident that the continuum slope did not change in the span of three weeks. Fewer emission lines were detected in the second spectrum, obtained with cloudy sky and thus noisier than the first one. The stellar absorption lines, if present, are embedded in the noise. The S/N value for the continuum at λ = 5500 Å is 180 and 45 in each spectrum, respectively. The EWs measured in both RSS spectra have approximately the same values; however, the EW of the hydrogen lines is smaller in the HRS spectra obtained in 2017 than in the 2016 RSS spectra, while the EWs of the He i lines show a marginal increase (see Table 2).

We found that the radial velocity is in the 150–190 km s⁻¹ range, which is consistent with SMC membership (systemic velocity varying from 88 to 215 km s⁻¹, according to Stanimirovic et al. 1999). In the RSS spectra, we measured larger values of Δv than in the HRS spectra. If we compare the σ and Δv values in Table 3 with the EW values in Table 2, we infer a slight increase of disk velocity from both Hβ and Hα between 2016 and 2017 corresponding to a decrease in EW values, as expected (Hanuschik et al. 1988; Reig et al. 2016), while there are no significant variations within timescales of weeks.

Also, for SUZAKU J0105–72, we were able to confirm the presence of a Be star, although the spectra did not allow detection of many emission lines. We detected Be-typical Balmer hydrogen and helium emission and absorption lines on two different dates with a 29 day span in between. The weather conditions were cloudy during the first run, so the resulting spectrum is weaker and noisier than the second one (see Figure 2, with S/Ns of about 25 and 60, respectively, measured at λ = 5500 Å). The two low-resolution spectra show similar continuum slopes, and the small variations in the red and blue portions of the spectral range (of less than 10%) are due to small differences in the sensitivity curves used for the photometric calibration. In both spectra, Hβ is not clearly detectable because it is stronger in absorption than in emission. The absorption should also affect the Hα emission line; therefore, our EW measurements for this object are lower limits. The Hα EW in the first low-resolution spectrum is lower than in the second, while high-resolution spectra show similar values. We measured radial velocity in the range 160–190 km s⁻¹ in the RSS and HRS observations, which is consistent with 186 km s⁻¹ measured in a catalog by Evans & Howarth (2008). We observed a small variation of Δv in 2017 between July/August and October for Hβ. In 2016, the σ values of Hα did not seem to have varied. In 2017, we found variations between August and October, while the EW values of Hα and Hβ are approximately constant. Therefore, for this object, the expected relation between EW and Δv does not seem to be valid. Although small variations of EW are expected in the case of small variations of Δv, according to the Hanuschik et al. (1988) relation, we do not have sufficient measurements to draw a definite conclusion. Reig et al. (2016) used data obtained over 15 yr in their Figure 7, so we suggest that much longer and possibly denser monitoring is necessary to draw conclusions.

XMMU J052016–692505 was observed only once with the RSS and the PG0900 grating. Our spectrum is of better quality than the one of Kahabka et al. (2006), and, in addition to the emission lines of hydrogen, we detected He i and Fe ii emission lines, making the Be star identification more solid. We also found that in this spectrum, the hydrogen emission lines are affected by absorption so strong that Hδ and Hγ are not measurable; therefore, our measurements represent only lower limits for the EW of Hα and Hβ. The source was also observed twice with a span of one month in between with the PG2300 grating.

The Hβ profiles are similar and show a red asymmetry. We fitted them with two Gaussian functions. In both spectra, the Hβ emission line is embedded in the corresponding absorption line. The radial velocity (v ∼ 280–300 km s⁻¹) is consistent with LMC membership. For Hβ, we detected a small decrease of Δv values between 2016 and 2017 and a small increase for the Hα line. The EW and Δv of Hα seem to follow the Hanuschik et al. (1988) relation, even if we have only four data points. For Hβ, we have more data points, but the correlation is less definite, even if there is an indication of decreasing EW for higher values of Δv.

Because the published spectrum of MAXI J0158–744 (Li et al. 2012) is of good quality, we focused on improving the resolution on the spectral range of Hβ and the He ii line at 4686 Å using only the PG2300 grating. It was observed twice with a 13 day span in between (Figure 3). The spectral range is centered on the Hβ emission line, which is asymmetric with a blue wing in both spectra. The Fe ii λ5018 is detected in our RSS spectra, but it was not measured by Li et al. (2012). It is noisy and asymmetric in both of our spectra. The weak He ii emission line at 4686 Å, detected by Li et al. (2012), is not observed in our spectra. The range of radial velocity is between 140 and 200 km s⁻¹. The Hβ profiles do not vary much, and we fitted them with two Gaussian functions. For both Hα and Hβ lines, Δv does not show variations. The EWs do not vary for both Hβ and the Fe ii emission lines in the RSS spectra, but larger values were measured for Hβ in the HRS ones. The Hα and He i EW values are consistent within the errors. The EW and Δv of both Hα and Hβ do not follow the expected relation, but we also do not have enough data to draw a definite conclusion.

4.2. Analysis of the Emission-line Profiles

The higher-resolution (PG2300 and HRS) spectra allow us to study the emission-line profiles and compare them at different epochs. We plotted the profiles in velocity space in Figures 4–7, assuming as the systemic velocity of the SMC sources the average radial velocity of the SMC, 158 km s⁻¹, and we adopted 278 km s⁻¹ for XMMU J052016.0–692505, which is in the LMC.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

XMMU J010147.5–715550 was observed with the HRS at three different epochs. The Hβ line shows two peaks with V/R ∼ 1 in the first-epoch spectrum, then the ratio increased and decreased again. The Hα line does not show a distinct double-peaked profile, and it is quite different from the Hβ. It is asymmetric, with three peaks in the first epoch, but with a "wine-bottle" profile with two peaks in the second and third epochs (see Figure 4).

The SUZAKU J0105–72 HRS spectra are noisier than those of the other sources but show weaker hydrogen emission lines than in the spectra of other sources. Nonetheless, both Hβ and Hα have a clear double-peaked profile, which is almost flat at the first epoch but still allows a double Gaussian fitting, at least for Hβ. The V/R ratio for Hβ decreased between 2017 July and two observations in 2017 October, and in 2017 October, we could also measure it for Hα, finding agreement with the Hβ values. We notice that the profiles of both lines underwent a similar evolution through the three epochs (see Figure 5).

For XMMU J052016–692505, we plotted the Hβ line observed with PG2300 in 2016 and both Hβ and Hα for the HRS 2017 observations. This is the most interesting object of our sample because the emission-line profiles show a clear evolution (see Figure 6). In 2016, Hβ had an asymmetric profile with a red wing, possibly due to a less bright, unresolved R peak. In 2017, the profile had changed shape with a red peak and a blue wing, which is in fact a blue peak of diminished intensity, rather stable for the rest of 2017. The V/R values decreased between 2017 July and the two 2017 October observations, and the Hα showed evolution in 2017. While two peaks may be present but are difficult to resolve in the first epoch, a triple-peaked profile appears in the second and third epochs. The general shape of Hα resembles that of Hβ, with a less intense blue peak, and it is possible that we did not resolve three peaks for the Hβ line.

Also, for MAXI J0158–744, we took both PG2300 and HRS spectra for Hβ. There is little evolution of the line profiles, which show a red peak permanently brighter than the blue one. Also in this system, the Hα line profiles are different from the Hβ profiles. The profiles are all asymmetric with a clear excess in the blue side. Three peaks can be distinguished in the first Hα spectrum, although only two are prominent, and the same may be true in the last Hβ spectrum.

We have confirmed the Be classification of all four optical counterparts of these SSSs whose spectra have the typically double-peaked emission lines of Be stars (see, among others, Struve 1945; Dachs et al. 1992; Hanuschik 1996; Labadie-Bartz et al. 2018) and radial velocities that are consistent with Magellanic Cloud membership. These four objects are clearly similar to the "emission-line Be stars" classified by Porter & Rivinius (2003), rather than to the "shell Be stars," because the central reversal point between the peaks of the Balmer hydrogen lines and some other lines of He i and Fe ii is always above the flux level of the stellar photosphere (however, the shell Be may only be an effect of inclination; see Hanuschik 1996). It is interesting that the excretion disk seems persistent over the years; in many Be stars, the emission-line spectrum is not observed at all epochs. Thus, the hypothesis of Sturm et al. (2012) of a transient excretion disk for XMMU J010147.5–715550 seems unfounded.

We have not detected any optical emission line indicating high ionization or excitation potential; most notably, the He ii λ4686 line is missing, showing that, most likely, there is no accretion disk. Not many clear-cut cases of accretion disks in Be binaries are known, but one is a rare Be+black hole system, MWC 656, and in that system, the He ii λ4686 line is clearly detected (Casares et al. 2012, 2014). Another candidate Be+black hole system is ϕ Per (Poeckert 1979), but this is not even a detectable X-ray source. However, the perfect spatial coincidence of five Be stars with SSSs cannot be coincidental, so these Be stars should have a binary compact object companion from which the SSS emission originates. We suggest that the optical spectra imply that these binaries are similar to the many neutron star X-ray binaries in which evidence of an accretion disk is missing: the neutron star accretes either when it crosses the excretion disk of the Be star or when there is renewed activity of the Be that increases the disk size and amount of material. In several known high-mass X-ray binaries (HMXBs), there is evidence that a temporary accretion disk is formed that does not remain for long (Reig 2011), and although we cannot rule out that the accretion disk is ever formed for the systems studied here, it seems very unlikely that a disk feeding a black hole in a super-Eddington accretion regime can have been formed and disappeared in such a short time that it was never observed for our sources. There are some interesting observable properties in these objects that should be monitored in the future. We discuss here the points that are relevant for future programs.

(1) In the absence of an accretion disk, the velocity difference Δv between the red and blue peaks of the most prominent emission lines that we obtained from the Hβ line, and in one case also for the Hα line, is not indicative of orbital motion; it indicates the excretion disk rotation velocity, as commonly assumed for Be stars. Although we measured variations of the line profiles and V/R ratio within weeks, Δv varied within a few weeks only for SUZAKU J0105–72. In the other objects, velocity variations occurred on a 1 yr timescale. The derived velocity is much lower (Δv = 100–200 and 100–220 km s⁻¹ for Hα and Hβ, respectively, corresponding to v_d sin i of about 50–100 and 50–110 km s⁻¹) than the average rotation velocity of galactic Be stars of different types (223 km s⁻¹ for early-Be type with weak emission lines, 270 km s⁻¹ for all others; see Briot 1986). Hanuschik (1996) found that in many Be stars, the disk is differentially rotating, with its axis aligned with the stellar rotation axis and a velocity law that is radially decreasing. Indeed, we know that the disks of both classic and binary Be stars are in Keplerian rotation (Meilland et al. 2012; Rivinius et al. 2013), but the latter are denser than the former and truncated (Reig et al. 2016). In a disk in Keplerian rotation, the radius is proportional to the inverse square of the velocity. If R_* is the radius of the Be star and v_* is its rotation velocity, the disk radius r_d can be inferred from Δv as

$$\begin{eqnarray*}&&({r}_{{\rm{d}}}/{R}_{* })={(2{v}_{* }\sin i/{\rm{\Delta }}v)}^{1/j},\end{eqnarray*}$$

where, of course, 1/j = 2 for Keplerian rotation (Reig et al. 2016). By applying this formula and using our mean value Δv = 150 km s⁻¹ and the mean value of v_* sin i from Briot (1986; 250 km s⁻¹), we obtain r_d ∼ 11 R_*, implying a relatively small disk compared to the values for classic Be stars (14–22 R_*; Reig et al. 2016).

(2) We found that the emission-line profiles change in time, but while for SUZAKU J0105–72, Hβ and Hα have similar profiles in the same spectrum, in two other sources, XMMU J010147.5–715550 and MAXI J0158–744, each of these two emission lines is different from the other at a given epoch. In MAXI J0158–744, there are only small differences in the line profiles at a given epoch. We note that in XMMU J052016–692505, there is even an inversion in the height of the blue and red peaks in Hβ. Because Hα has a lower excitation potential, it is supposed to arise in a cooler region than Hβ. According to Sigut & Jones (2007), the thermal structure of the Be disk involves a higher temperature in the "edge" region at a high azimuthal height. Moreover, the optical depth of Hα is higher than that of Hβ, so the Hα line appears stronger (forming over a larger radial portion of the disk). In a recent article, Panoglou et al. (2018) calculated models of excretion disks of binary Be stars with radiative transfer and found that binary tidal effects cause phase-locked variability of the V/R ratio of the Balmer lines, with two maxima per orbital cycle. Thus, we may be able to infer the orbital period of our systems by measuring the V/R with a short cadence of one or two weeks, perhaps.

(3) Reig (2011) and Reig et al. (2016) found interesting correlations of the EW of the Hα line in all Be stars with excretion disk size and HMXBs with the orbital period. Reig (2011) proposed that the orbital period determines the size of the disk, which is truncated by the tidal interaction of the secondary. Since the Hα line originates in the cooler outer zone of the disk, it is easily conceivable that its EW is inversely proportional to the disk size. Figure 15 of Reig (2011) shows that the maximum EW(Hα) in HMXBs is proportional to the orbital period. In our cases, for SUZAKU J0105–72, the orbital period should be of the order of less than 10 days (EW ∼ 4 Å), while it should be two to four months for our other three sources (EW in the range 22–30 Å). These are actually lower limits; in fact, the model assumes that the disk is dense enough to reach the tidal radius of the binary. This may not be true if there are periods of interrupted or decreased mass outflow from the Be star. In any case, even having upper limits on the period is useful as a benchmark to plan observations with a suitable cadence and measure the orbital periods of these binaries. However, in short orbital period systems, the EW may be variable. Reig et al. (2016) also studied the variability of the Hα EW in HMXBs and found that systems with short orbital periods are affected by the neutron star tidal truncation of the disk more often and more strongly; thus, the disk does not have a stable configuration, and its density changes quickly.

(4) Finally, we should also take into account two suggestions by Raguzova (2001). The first concerns the He i emission lines. These lines require larger ionization potential than the Balmer hydrogen lines, are typical of Be stars with a WD companion because they are formed in the Strömgren sphere near the WD, and indicate the WD motion. Raguzova (2001) calculated that, if the orbital period exceeds about three months, the He i lines should have a semi-amplitude of the velocity curve due to orbital motion of 70 km s⁻¹ that should be well measurable. Both this type of measurement and that of the broadening of the Be absorption lines to determine the rotational velocity of the Be star possibly require a larger telescope collecting area and certainly require longer exposures, more frequent sampling, and/or higher spectral resolution than those we used in the present project. Raguzova also predicted that the semi-amplitude of the velocity curve of the Balmer lines, associated with the excretion disk around the Be star, would only be 10–20 km s⁻¹. Of course, in the Magellanic Clouds, the measured radial velocity of the Balmer lines is due to the systemic velocity in the region of the Clouds in which each star is located plus the component due to the orbital motion. To obtain the semi-amplitude from the velocity curve, higher spectral resolution and measurements repeated with more frequency (every week or once a day) are needed.

In the previous section, we outlined possible photometric and spectroscopic observations that may yield measurements of the orbital period and WD rotation period and would allow us to derive the masses of the Be star and the compact object. While certainly there is more interesting work ahead before completely understanding the nature of these systems, we are not aware of a phenomenology that would make neutron stars intermittent or transient supersoft X-ray emitters without also presenting harder X-ray states and/or a much harder X-ray "tail" in addition to the SSS spectrum (see Kylafis & Xilouris 1993; Kylafis 1996; Kohno et al. 2000). The SSS phenomenology seems rather typical of burning WDs or, more rarely, of black holes with extended accretion disks. However, the Be+black hole scenario is also highly unlikely: large, luminous, steady, and "cool" disks, like the one in the ultraluminous SSS in M101 (Liu et al. 2013), and perhaps in the M81 source (Liu et al. 2015), do seem to produce a supersoft X-ray spectrum, but they are persistent SSSs, and there is a prominent He ii line in the case of M101. We suggest that the absence of this emission line and the short duration of the SSS phases in our systems cannot be reconciled with an accreting black hole as compact object.

The clear association of five Be counterparts out of about 30 SSSs in the Magellanic Clouds suggests that such binaries may be common. The fact that many SSSs in M31 and galaxies outside the Local Group are in star-forming regions (Di Stefano & Kong 2003; Orio et al. 2010) also indicates that some transient or recurrent SSSs may be the detectable manifestation of Be stars with WD companions, thanks to the luminous X-ray property of shell hydrogen burning, ignited perhaps intermittently, as the X-ray observations of these four sources seem to indicate. Though difficult to discover, Be+WD star systems may be common according to the works we mentioned in the Introduction, so this may be a whole new class of interacting binary progenitors of SNe Ia.

The model of massive WDs accreting from Be stars and burning hydrogen has been proposed to explain several cases of SSSs, but the evolution of such systems has not been explored yet. Would they eject the accreted material in a nova-type outburst? A nova in a Be star system may easily be missed, having an amplitude of less than 3 mag, not even comparable with the optical novae we know in low-mass systems. Moreover, if the WD accretes and burns only sporadically during the orbital cycle, mass ejection may be rare or impossible, and the already massive WD could grow toward the Chandrasekhar mass. There is one more parameter to take into account, and it is the WD temperature at the beginning of mass accretion. Because the WDs in these systems are very hot, they have less degenerate accreted envelopes, and most likely, even if they undergo nova outbursts, they eject less mass in nova events and do not undergo phases of helium shell flashes (Hillman et al. 2016), possibly continuing to accrete more material than they can ever eject and being on a path toward SN Ia explosion– or accretion-induced collapse into neutron stars.

The dependence of the SN Ia rate on star formation rate has been known since the pioneering works of Tammann (1977, 1978) and Oemler & Tinsley (1979) and has been confirmed and explored time and again by different authors (see Mannucci et al. 2006; Sullivan et al. 2006); it will be very interesting to explore how this new class of progenitors fits into galactic evolution and explains SNe Ia in young populations. In evaluating the statistics and the time for accretion toward the Chandrasekhar mass, we must remember that the Be+WD binaries are not long-lived, but they last longer than the numerous neutron star+Be star systems we know, and their short-lived SSS phases make them difficult to identify even in X-rays.

Following the work of Panoglou et al. (2018), the changing profile of the optical emission lines in Be stars is an indication of binarity. We also note that the appearance of a triple-peaked profile, observed sporadically in our spectra, is not a proof against binarity, according to the above authors. In fact, the Be binary zeta τ does present triple peaks (Carciofi et al. 2009). Perhaps the way to detect the presence of a WD in a Be star without a known orbital companion may be cyclic variability of V/R (due to the orbital period of the otherwise undetectable WD). This would be a more promising way of discovering Be+WD systems than "waiting" for an SSS flare. If we want to build statistics of Be+WD binaries, the X-rays probably only show the tip of the iceberg, while optical studies are a more feasible avenue to increase the number of known systems.

All of the observations reported in this paper were obtained with the Southern African Large Telescope (SALT). We are grateful to the anonymous referee for useful comments and suggestions that helped to improve the quality of the paper.