© ESO 2022

1. Introduction

The proposed stellar mass black hole (BH) candidates in the long period binary LB-1 (Liu et al. 2019) and the triple system HR 6819 (Rivinius et al. 2020) are possible additions to the class of X-ray quiet Be+BH systems of which MWC 656 (Casares et al. 2014) is the only confirmed member to date. These systems, with periods of some tens of days, are especially interesting as theory (Langer et al. 2020a) predicts that they may contribute to the population of BH+BH systems that led to the discovery of gravitational waves (Abbott et al. 2016). However, the nature of these new systems is disputed (Abdul-Masih et al. 2020; El-Badry & Quataert 2021; Simón-Díaz et al. 2020; Irrgang et al. 2020; Shenar et al. 2020a; Bodensteiner et al. 2020; Lennon et al. 2021).

Morphologically, the spectra of LB-1 and HR 6819 are composite, having a narrow-lined B-type spectrum superimposed on a relatively featureless Be-type spectrum from which only the hydrogen Balmer and stronger He I lines are detected. The B-type spectrum displays clear radial velocity variations of some tens of km s⁻¹ consistent with a long-period binary, while the Be-type spectrum is almost stationary to within a few km s⁻¹. Interpreting the latter as binary motion anticorrelated with that of the B-type star leads to mass ratios on the order of 10, and the systems’ natures then hinge on the derived masses of one or another of the components. If the B-type star is a main sequence star, then the Be-type spectrum may signify a candidate BH or neutron star (NS) plus accretion disk (Liu et al. 2019; Abdul-Masih et al. 2020), though tension may arise if the mass of the B-type star is small enough such that the implied reflex motion of the BH or NS accretion disk becomes larger than observed (Simón-Díaz et al. 2020). Conversely, if the Be-star is a classical Be star, then the B-type star may be a stripped helium star contracting to become a hot subdwarf (Irrgang et al. 2020; Shenar et al. 2020a; El-Badry & Quataert 2021; Bodensteiner et al. 2020). Finally, the Be-like emission can be interpreted as being stationary (El-Badry & Quataert 2020; Abdul-Masih et al. 2020) and unrelated to the secondary of the B-type star, either belonging to a third very long period Be-star in the system (Rivinius et al. 2020), or even as chance alignment with a Be-star (Safarzadeh et al. 2020). In this case, the secondary to the B-type star is not identified in the spectrum. In the context of triple systems hosting BHs, HD 96670 may also be relevant as Gomez & Grindlay (2021) suggest it may be a triple system with an inner short-period (∼5.3 d) binary consisting of a 23 M_⊙ O-star plus a 6 M_⊙ BH.

The star considered here, with the Simbad designation Cl* NGC 2004 ELS 115 (hereafter NGC 2004#115), lies on the periphery of the young cluster NGC 2004, on the northern edge of the Large Magellanic Cloud (LMC). It was observed as part of the VLT-FLAMES Survey of Massive Stars (FSMS, Evans et al. 2005, 2006) and classified as a B2e star (without a luminosity classification). It was noted to be a single-lined spectroscopic binary (SB1) with double-peaked Hα emission that also has a central absorption component. The star was subsequently discussed by Dunstall et al. (2011) in their study of nitrogen abundances in Be stars from the FSMS, where they noted it as a short-period system and likely not a classical Be star. They gave its projected rotational velocity (v_e sin i) as 15 km s⁻¹, which is close to the spectral resolution of the data, but it was not considered further.

However, NGC 2004#115 recently attracted our attention as the morphology of its spectrum is strikingly similar to that of LB-1 and HR 6819. It has a narrow-lined B-type spectrum that is clearly undergoing short-period (a few days) radial velocity variations. There is no clear evidence for a second spectrum other than the Hα emission line that is almost stationary but it may exhibit small antiphase variations in the two epochs obtained (see Fig. 1 as well as Sect. 3.3 and Fig. 10 for a further discussion of Hα). However, additional observations were recently obtained with the Southern African Large Telescope (SALT; Buckley et al. 2006) to improve phase coverage and revealed the Hα emission to be much weaker, and provide strong evidence that NGC 2004#115 is a triple system.

Fig. 1. Composite nature of the Hα profile that consists of a broad double-peaked emission feature, plus a narrow absorption component from the B-type star, the narrow lines of the C II 6578/6582 Å doublet are further to the red. The two epochs from the VLT-FLAMES observations shown, OB1 (red) and OB2 (blue), illustrate significant velocity variability of the narrow components but a relatively stable broad Hα emission.

2. Observational data

The FSMS data are presented in detail by Evans et al. (2006), therefore we only summarize the relevant points here. NGC 2004#115 was observed with VLT-FLAMES between late 2003 and early 2004 with each of six high resolution grating settings; HR02 (central wavelength 3958 Å), HR03 (4124 Å), HR04 (4297 Å), HR05 (4471 Å), HR06 (4656 Å) and HR14 (6515 Å)¹, with wavelength coverage of each region listed in Table A.1. All regions were observed in two observing blocks (OB1 and OB2), each of which was split into 3 contiguous 2275 s exposures, for a total of 36 spectra with a typical resolution of ∼20 000². However the signal-to-noise ratio (S/N) for the second OB of the HR03 setting was poor for all stars in this configuration, with only the first observation proving to be useful.

NGC 2004#115 was also observed with SALT using the High Resolution Spectrograph (HRS; Bramall et al. 2010). The four observations (17, 19, 22, and 23 September 2020) were performed in low resolution mode (R ∼ 14 000) with exposure times of 2400 s and covered a wavelength range of approximately 4000–8800 Å. Primary reductions were executed with the SALT science pipeline (Crawford 2015) which include overscan correction, bias subtraction and gain correction. The subsequent reduction steps, which include background subtraction, arc line identification, removal of the blaze function, and merging of orders, were carried out with the MIDAS FEROS and ECHELLE packages. Orders containing the broad Hα and Hβ lines were reduced following the recipe discussed in Appendix A of Simón-Díaz et al. (2020).

3. System properties

The data for this system exhibit the spectral signatures of three components; a narrow lined B-type spectrum of spectral type B2 (illustrated in Fig. A.1), an emission-line Be disk-like spectrum characterized by Balmer emission that is variable (no other emissions lines are detected at the present S/N), and the faint signature of another star in wings of the Balmer and He I lines, as shown in Fig. 2. We tentatively identify this last component as belonging to a third star, and suggest that the close companion to the narrow-lined star is undetected. For convenience we refer to the narrow-lined B-type stars as the primary (sometimes as m₁), its undetected close companion as the secondary (m₂), and the faintly detected star as the tertiary (m₃). The source of the emission line spectrum we refer to as the disk. We refer to the primary+secondary system as the inner binary (at times the binary) and the inner binary plus tertiary as the outer binary.

Fig. 2. Left panel: comparison of Hγ profiles from two VLT-FLAMES epochs in the rest frame of the narrow-lined B-type star. At these epochs the star is close to each quadrature phase. An arrow indicates the extra absorption that we attribute to the presence of a third star, not present in the blue spectrum. If this component was due to a faint companion to the B-type star then it should be visible on the red wing of the blue spectrum. Instead its velocity at that epoch is almost identical to that of the B-type star, see Sect. 3.2, hence no absorption is apparent at that epoch. Right panel: as for the left panel except here we plot two VLT-FLAMES epochs of the He I 4471 Å line. Similar considerations apply.

3.1. Properties of the narrow-lined B-type star

3.1.1. Radial velocity analysis

The sharpness of the metal lines in the primary enables precise measurement of its radial velocity, described in Appendix A, with the results listed in Table A.2. The VLT-FLAMES velocities were initially used as inputs to the adaptive simulated annealing code rvfit (Iglesias-Marzoa et al. 2015) to solve for the orbital parameters of the primary. This resulted in the determination of orbital parameters with period (P) of 2.88 days, eccentricity (e) of 0.27, and systemic (γ) and semi-amplitude (K₁) velocities of 301 and 74 km s⁻¹ respectively. The residuals were small, less than 6 km s⁻¹, although the 1-σ errors on the radial velocities were estimated as ∼1 km s⁻¹. However addition of the SALT data resulted in a similar solution but with worse residuals, up to ∼10 km s⁻¹, and with the overall rms of the fit increasing from 3.5 to 4.6 km s⁻¹. This prompted a more careful examination of the VLT-FLAMES data. Noting that approximately two thirds of these data were taken within a one week period a solution was sought with these data only. This resulted in a much better fit to the data, though with negligibly small eccentricity. The RVFIT orbital solutions for these three cases are shown in Fig. 3, while the detailed results for this last solution are listed in Table 1.

Fig. 3. Radial velocity (RV) orbital solution for NGC 2004#115 based on velocity estimates from the metal lines in its spectrum. Left-most panel: illustrates the fit using the VLT-FLAMES data only. Central panel: how the fit, now determined using both VLT-FLAMES (black points) and SALT-HRS (red points) data, deteriorates when we include both datasets. Rightmost panel: fit to the VLT-FLAMES data taken within the final week of that observing campaign, and fit parameters for this solution are listed in Table 1. The scales differ from panel to panel, in particular note the change in the scale of the O−C plots.

Inspection of the fit to the early VLT-FLAMES and SALT-HRS data hints at a modulation of the radial velocity curve on longer timescales, though with small amplitude, as illustrated in Fig. 4. An obvious explanation is that NGC 2004#115 is a triple system, and indeed there is additional evidence from the SALT-HRS data of a blended component to the Hβ line (discussed below). From the radial velocity discrepancies the semi-amplitude velocity for this tertiary component appears to be at least ∼20 km s⁻¹, while the trend of the VLT-FLAMES data (early and late) over a 60 day time-span has no obvious minimum or maximum arguing for a period in excess of 120 days. However, given the sparseness of the sampling of the outer binary these are extremely tentative statements.

Fig. 4. Central panel: comparison of the solution derived from the late epoch FSMS data only (solid line) with the observed radial velocities of the primary, filled circles plus (small) error bars. Left panel: comparison of that solution with the early epoch FSMS data. Right panel; comparison with the SALT/HRS data. Horizontal lines indicate the systemic velocity of this solution. The dotted lines in the right panel are the solutions obtained by fitting those radial velocities but fixing both the period, eccentricity and semi-amplitude velocity to the values derived from the late FSMS data. Open star symbols are the velocities of the emission features (Table 5), and open circles are velocities of the broad-lined component (Table 4) that we identify as a possible tertiary.

3.1.2. Atmospheric analysis

The grid of non-local thermodynamical equilibrium (NLTE) models used in this analysis was calculated with the TLUSTY and SYNSPEC codes (Hubeny 1988; Hubeny & Lanz 1995; Lanz & Hubeny 2007) and adopted the “classical” model atmosphere assumptions, that is plane-parallel geometry, hydrostatic equilibrium, and that the optical spectrum is unaffected by winds. They also assumed a normal helium to hydrogen ratio (0.1 by number of atoms) and element abundances appropriate to the LMC. Microturbulence was not included as a pressure term in the hydro-static equilibrium equation. Further details can be found in Ryans et al. (2003) and Dufton et al. (2005, 2018), while the analysis methods are discussed further in Appendix B.

As already noted, there is little evidence of a second spectrum present other than the Balmer emission, and potentially broad blended Balmer lines, both of which are discussed below. However, the sharp metal lines are well observed except for the veiling effect of other sources. Indeed a straightforward analysis of the spectrum assuming no veiling leads to a solution in which all elements are underabundant (discussed below and illustrated in Table 2). Therefore, following Lennon et al. (2021), for a given log g the degree of veiling can be estimated by requiring that an abundance analysis results in a normal silicon abundance for the LMC. We adopt the LMC silicon abundance of 7.2 ± 0.1 dex from Hunter et al. (2009) as they methods similar to those used here, thereby minimizing systematic effects. In fact many of their sample of B-type stars are also in the NGC 2004 cluster, and its vicinity.

The stellar effective temperature (T_eff) is relatively well constrained by the presence of Si III lines, the near absence of Si II and Si IV, and the absence of He II lines. As shown in Fig. 5, the contribution of the narrow-lined star to the total flux in the V-band lies in the range 0.54–0.62, with T_eff  ∼  21 500 ± 1500 K depending on log g, and assuming a helium abundance by number of N[He/H] = 0.1. Increasing the helium abundance to N[He/H] = 0.2 decreases the flux contribution from the B-type star by ∼0.02–0.04 but, as we find that the range of surface compositions implies only a slight nitrogen enhancement, a helium enrichment is unlikely. Therefore we have not explored helium-rich models further. The uncertainty in [Si/H] of ±0.1 translates into an uncertainty of ∓0.05 in the flux contribution from the B-type star.

Fig. 5. Locus of (log g, Teff) parameters, labeled with the fractional B-star contributions, that reproduce the normal LMC silicon abundance of 7.2 dex. A helium abundance of N[He/H] = 0.1 and a microturbulence of vt = 0 km s−1 have been assumed.

For early type stars, the primary diagnostic of the surface gravity is the strength of pressure broadened wings of the hydrogen Balmer lines. In the present case, the composite nature of spectrum prevents use of this standard method, while a disentangling approach (as in Shenar et al. 2020b) is also not feasible since, for example, there are effectively only two epochs of the Hγ observations. However, the Balmer lines, or the gravity sensitive Diffuse He I lines, offer an alternative method for determining the degree of veiling, and by inference log g, as their sharp Doppler cores provide a method for assessing which of the models indicated by Fig. 5 best removes that feature from the residual spectrum. The Balmer lines are not used since modeling the line cores themselves depends on the mass-loss rate (Najarro et al. 1996), which is not included in our models. By contrast, the weaker He I lines are well modeled by NLTE plane parallel hydro-static equilibrium codes (see for example Irrgang et al. 2020).

Therefore, as the helium abundance is essentially solar (from above), we use the He I 4471 Å line to constrain log g. This is the strongest Diffuse line in our data, it has the best understood line broadening, and lies in the region of highest S/N of the VLT-FLAMES data. Furthermore, from examination of FSMS data for other stars within ∼30″ of NGC 2004#115, ELS numbers 23, 73, 98 and 118, it is also clear that nebular emission is insignificant. For this test we merged the subexposures within each observing block of the VLT-FLAMES data, with appropriate velocity corrections. The results, indicated in Fig. 7, demonstrate that the residuals of the subtracted Doppler core imply a log g ≥ 3.75, though the upper bound is poorly constrained. We therefore adopt log g ≃ 4.0, and hence T_eff = 22 600 K and a B-type star flux contribution of 0.58. We recognize that log g is uncertain, however, as discussed in Sect. 4, the star’s luminosity, as derived below, provides a more precise evolutionary mass, which is consistent with log g ≃ 4.0 (the determination of the radius and luminosity is insensitive to the value of log g).

Using the above information, and comparing the observed de-reddened V-band flux (Sect. 3.4) with that computed using the filter function of Bessell (1990), and assuming a distance to the LMC of 49.59 kpc (Pietrzyński et al. 2019), leads to estimates of the stellar radius and luminosity of the B-type star as indicated in Table 2. This table also indicates the derived chemical composition and, for comparison, equivalent results for the case of zero veiling. The composition of the star is very similar to baseline LMC abundances, with perhaps a modest nitrogen enhancement being apparent. As discussed in Appendix B, while we adopt a best fit microturbulence velocity (v_t) of v_t = 0 km s⁻¹, increasing this to v_t = 2 km s⁻¹ would strengthen the theoretical metal lines and lead to a decrease in the B-star contribution to ∼0.51, for the above value of log g.

3.2. The properties of the broad-lined star

Adopting the previously determined parameters for the primary of T_eff = 22 600 K, 0.58 fractional contribution to the V-band flux, helium abundance N[He/H] = 0.1, and log g = 4.0 dex, we subtract the model spectrum from the data, appropriately weighted and Doppler shifted, to produce a residual spectrum. To improve the signal-to-noise we merge individual exposures per OB of the VLT-FLAMES data. The resultant spectra exhibit He I lines that support our assumption that the tertiary is a mid to early B-type star, see Fig. 6, moreover we can estimate its radial velocity from the He I and Balmer lines (noting that the presence of emission may impact the results).

Fig. 6. Subtracted spectra of the proposed tertiary for the HR02, HR04 and HR05 regions. We have not included HR03 due to its low signal-to-noise (discussed in the text), or HR06 as it is relatively featureless. Close examination of this spectrum may show reveal the presence of over-subtracted metal lines, but these are artifacts due to the use of a specific model spectrum from our grid with mean surface abundances a little different from the final abundance determination of the primary. Spectra were approximately normalized using a low order polynomial fitted to continuum regions.

Fig. 7. Illustration of the residual spectra obtained by subtracting theoretical spectra scaled by the fractional contribution from the HR06 OB2 data in the vicinity of He I 4471 Å line. The value of log g increases in steps if 0.25 dex from log g = 3.00 dex for the lower curve to 4.25 dex for the uppermost curve, spectra are offset by 0.2 normalized flux units. Teff values follow the relationship shown in Fig. 5. The emission feature at the position of the He I line core (dashed line), caused by over-subtraction of the Doppler core in the model, becomes weak or absent for log g > 3.75.

In determining radial velocities and v_e sin i estimates we used the He I lines at 4009/4026 Å (both in HR02 spectra) and 4387/4471 Å (both in HR05 spectra), other lines were found to be too weak in the residual spectrum to provide useful measurements or a confident detection. We followed the methodology of Dufton et al. (2022), who determined radial velocities, and v_e sin i estimates, for a large sample of Be stars in the LMC and SMC using the He I lines. Briefly, v_e sin i was estimated for each line using two different methodologies, profile fitting (PF) and Fourier transform (FT), and are listed in Table 3. The overall means for all lines agree quite well, but there are some differences between the PF and FT estimates per line and between lines; we adopt v_e sin i ∼ 300 km s⁻¹. Radial velocities were determined from the centroids of the lines resulting from the PF approach. The laboratory wavelengths for the allowed components of the He I 4026 and 4471 Å lines were corrected for the presence of forbidden components in their blue wings, based on simulations for an extensive grid of rotationally broadened model spectra. The radial velocity estimates are listed in Table 4 and, while formal errors are not given, there are differences of up to 30 km s⁻¹ between lines in a given VLT-FLAMES observation. This probably reflects the difficulty in measuring line centroids of the He I lines in such a broad lined star, but may also be partly due to a weak emission component (presumably from the disk or faint diffuse nebulosity) that is distorting the line profiles in the VLT-FLAMES data.

We also measured the radial velocities of the residual Hγ profiles from the VLT-FLAMES data (two epochs of HR04 data), and the Hβ profiles from the SALT-HRS data (all four epochs), though for these we fit rotationally broadened model profiles to the measurements, also listed in Table 4. The Hβ line is particularly well modeled as we can explicitly account for the weak and narrow emission component in fitting the residual spectrum (next subsection, and Fig. 8). These results in particular imply an approximate mean velocity of +267 km s⁻¹, compared to the systemic velocity of the binary of +320 km s⁻¹ at this epoch, and show very little evidence of variability despite spanning six days and very different phases of the inner binary. Evidently the residual broad-lined spectrum is that of the tertiary, and is not a component of the short period inner binary. For convenience the broad absorption line velocities are also plotted in Fig. 4.

Fig. 8. Residual (i.e., primary subtracted) Hβ profiles from the SALT/HRS data, observations 1 through 4 are from bottom to top, shifted by 0.5 units for clarity. The vertical lines indicate the mean positions of the absorption components (dotted line) and emission components (dashed line), from Tables 4 and 5.

3.3. The emission line spectrum

As for our discussion of the residual spectrum above, the following analysis of the emission features utilized the merged spectra for each OB of the VLT-FLAMES data. Broad double-peaked Hα emission is observed in the VLT-FLAMES data, but in the SALT-HRS data it is much weaker, narrower, and rather flat-topped, with just a hint of a double peak (see Fig. 9). Also apparent is a very weak contribution from diffuse nebular emission, best seen in the VLT-FLAMES spectra in Fig. 1 as a weak and narrow spike in the core of the emission line (its velocity is +315 ± 5 km s⁻¹). Diffuse nebular emission is observed at this velocity in the sky spectra, observed with SALT-HRS at distances of 40″ and 63″ from the source. The strength of this emission depends on position, however, and simply subtracting the sky spectrum from the source results in over-subtracted nebular emission, and strongly distorted line Hα line cores. Hence, the data shown in Fig. 9 have not been sky-subtracted and still exhibit the typical night sky emission lines observed near Hα. The SALT-HRS data reveal that the Hβ line, not observed by VLT-FLAMES, also exhibits weak emission. Very weak Hγ emission is detected in the VLT-FLAMES data, although its weakness prevented measurement of reliable emission line profiles.

Fig. 9. Hα profile morphology in the VLT-FLAMES data (top two profiles) and SALT-HRS data (four lower profiles). Spectra have been off-set by 0.5 normalized flux units for clarity. Sharp emission features in the SALT-HRS data are night sky lines (⊗). Hα and the sky line just long-ward of 6590 Å are likely blended with weak nebular H I and [N II] respectively, as discussed in the text. The position of the central dip of the Hα Doppler core of the narrow-lined B-type star is indicated, as are the positions of its two C II lines (↑).

There appear to be small velocity shifts in the broad Hα emission, as shown in Fig. 1, that are anti-correlated with those of the primary. However, as this apparent motion may be due to the neglect of the underlying absorption from the B-type star that distorts the emission profile (Abdul-Masih et al. 2020; El-Badry & Quataert 2020), we subtracted this contribution from the data, shown in Fig. 10. A small shift between the profiles is still observed and, given that to a good approximation they are symmetric outside of the central dip, we determine their radial velocities from the position of bisector of the profiles. This results in values of 310 ± 5 and 290 ± 5 km s⁻¹ for the merged OB1 and OB2 data respectively, see Table 5. The relative lack of motion of the emission is pronounced given the contemporaneous velocities of the B-type star of approximately +245 and +372 km s⁻¹ for OB1 and OB2 respectively and, while they appear to be anticorrelated, these measurements were taken 34 days apart and hence may reflect the systemic velocity drift of the inner binary.

Fig. 10. Left panel: VLT-FLAMES OB1 (blue) and OB2 (magenta) Hα profiles resulting from subtraction of the underlying spectrum of the narrow-lined B-type star. Horizontal lines, same color coding, represent the rest frame velocity of the primary at these epochs. Note that a small residual velocity shift is still apparent (cf Fig. 1). Right panel: the same profiles plotted in their rest frame velocities; OB1 and OB2 spectra in blue and magenta respectively, the merged spectrum is plotted in black. Also illustrated are the velocities of the double peaks (±155 km s−1), the full-width half-maximim (FWHM) of the emission (550 km s−1), and the extent to which emission is present in the line wings (±500 km s−1). These spectra were taken ∼34 days apart and some variability is apparent within the double-peaked structure of the profile.

The Hα emission in the SALT-HRS data is narrow, with just a hint of a double peaked structure (it as almost flat-topped) when corrected for the underlying absorption and night sky lines, see Fig. 11. Ambient nebular emission must also be present, for example a weak [O III] 5006.84 Å line is present in the SALT-HRS on-source data and we measure its velocity to be 314 ± 1 km s⁻¹, with a line width of 20 km s⁻¹ (the resolution of the data). Very weak [N II] lines are also present, with a consistent velocity. We also see weak Balmer emission with this velocity and width in the sky spectra, though as noted the sky fibers were quite far from the source position and the line strengths variable, preventing precise sky subtraction. The features shown in Fig. 11 are clearly broader than those of the ambient nebular emission, hinting at an additional emission component at higher velocity. We therefore fitted the Hα lines with a two-component Gaussian model, fixing position and width of one Gaussian at the values measured for the ambient emission, and optimizing the fit for its peak height and the parameters for the second component. The results for the velocity of this second component, which we identify with the source, are also summarized in Table 5, as are the Hβ results, for which we used a single component fit as the ambient feature there should be much weaker than in Hα. Clearly the unknown intensity of the ambient nebular emission is an issue, that could be resolved with high resolution spectroscopy. Nevertheless, the morphology of the emission is consistent with the source emission being at the higher velocity as the velocity of the ambient Hα emission is well constrained. From Table 5 we see that the Hα velocities, taken over a few days and at differing phases, are constant to within a few km s⁻¹. The weaker Hβ line, while it exhibits a hint of anticorrelation with the narrow lined star, has velocities consistent within the errors of the measurements. However, as already noted, the emission line velocities at this epoch are very different from the radial velocities of the broad-lined star (Table 4), and quite similar to the systemic velocity of the inner binary. On this evidence therefore we tentatively associate the source emission with the inner binary, and not with the broad-lined tertiary, although we discuss alternative scenarios in Sect. 5.

Fig. 11. SALT-HRS line profile morphology in velocity space of Hα (upper four curves) and Hβ (lower four curves), color coded by observation; black – obs1, blue – obs2, magenta – obs3 and orange – obs4. The primary’s underlying profiles have been subtracted from these data. While the Hα profiles exhibit some variability of intensity their velocities are quite constant, and are well modeled with two-component Gaussian fit with velocities at 314 km s−1 (for the ambient nebular emission) and ∼325 km s−1 (for the source emission). The latter has a FWHM of ∼60 km s−1.

3.4. Photometry, optical light curve, and X-ray flux limit

The available optical and infrared (IR) data are described in Appendix D, and in general are in good agreement with the results of the atmospheric analysis, though there is evidence of a small IR excess. The magnitude of the source has been extremely constant, illustrated by approximately 7 years of phase-folded MACHO data in Fig. 12, that limits possible ellipsoidal light curve variations to less than 0.005^m in V and 0.007^m in R in amplitude, or around 0.5–0.6%. The lack of variability of the system is also reflected in the Gaia data with the EDR3 G-band flux indicating an uncertainty of 0.1%, and simulating the modulation of ellipsoidal light variations in the data we estimate the amplitude of these variations to be less than 0.3%, see Appendix D for details. Finally, XMM-Newton observations provide an upper limit to the X-ray flux from the system of less than 5.33 × 10³³ ergs s⁻¹ (about 1.4 L_⊙) in the 0.2–12 keV band.

Fig. 12. MACHO V and R photometry spanning MJD 49001.6091 to 51542.5763, folded to the 2.92 d period of the system, illustrating the lack of variability in both bands. The magnitudes are instrumental; zero points and transformations to Kron-Cousins V and R magnitudes can be found in Alcock et al. (1999).

It the above context it is interesting that classical Be stars are well known variables, Hubert & Floquet (1998) find that 86% of early Be stars are subject to rapid variability of at least 0.02^m, equivalent to about 0.01^m given the veiling effect here, while NGC 2004#115 is constant to better than 0.002^m in the GaiaG band. In addition its IR excess is very small for a typical Be star, for example J − [3.6] = − 0.13^m compared with a typical value of −0.6^m (Bonanos et al. 2009, 2010).

4. Discussion

In this section we present some evolutionary scenarios and discuss to what extent they can explain the primary constraints on the system. We make the simplifying assumption that the orbital dynamics of the inner binary are not significantly perturbed by the tertiary, and adopt the mass function for this system from Table 1. Other hard constraints on this system are the ellipsoidal light variations (f_el) of less than 0.6–0.3%; the radius and luminosity of the primary from Table 2, its low v_e sin i < 15 km s⁻¹, and the lack of any eclipse that defines the maximum possible sini. The value of v_e sin i, together with the period and stellar radius, then define the sin i required to ensure synchronous rotation of the B-type star, or sin i < 0.15. Synchronous rotation would also imply that the equatorial rotational velocity of the B-type star is 97 km s⁻¹.

The mass of the primary (or m₁) is a critical parameter of the system that is not well constrained from observations (Sect. 3.1). The surface gravity of the primary implies a spectroscopic mass (M_sp) of ∼9 M_⊙, though with a large uncertainty giving a mass range of 5–16 M_⊙. However, the extremes of this mass range can be ruled out as being inconsistent with the star’s luminosity, for example, a 16 M_⊙ star would have a luminosity of 4.7–4.8 L_⊙, well above the estimate derived from its radius and T_eff. In fact we can also define an evolutionary mass (M_ev) using BONNSAI³ (Schneider et al. 2014, 2017) that is based on the models of Brott et al. (2011) and Köhler et al. (2015), which were calibrated using B-type stars observed by the FSMS, many of which are in NGC 2004 (Hunter et al. 2009). As input we use T_eff = 22 600 K, R₁ = 5.6 R_⊙, and v_e sin i as above, resulting in estimate of M_ev = 8.6 M_⊙ and an age of 22 Myrs. This is in reasonable agreement with the spectoscopic mass above, though better constrained, and consistent with the age of NGC 2004. Note that if this B-type star is a product of binary evolution then these mass and age estimates may not be applicable.

Alternatively the narrow lined B-type star might be a low mass stripped helium star (Wellstein et al. 2001; Götberg et al. 2018), with a more massive (∼10 M_⊙) Be-star companion, as has been proposed for LB-1 and HR6819 (Shenar et al. 2020a; Bodensteiner et al. 2020; El-Badry & Quataert 2021). In this picture, the stripped star is in a short-lived (∼10⁵ yr) transition as it evolves rapidly to its stable helium burning phase. Interpolating between the stripped star evolutionary tracks illustrated in the HR-diagram (HRD) in Fig. 3 of Langer et al. (2020b) we find that the luminosity of the narrow-lined B-type star of 3.87 L_⊙ corresponds to a stripped star mass of approximately 1.4 M_⊙, and with log g  ∼  3.2.

4.1. The narrow-lined B-type star as a stripped star

Considering first the stripped star scenario with m₁ = 1.4 M_⊙ we find that the Roche radius is smaller than the stellar radius by approximately 20%, irrespective of the value of sin i. The inflated stripped star is just too big to fit into such a short period system. The situation is illustrated in Fig. 13 which demonstrates that for primary masses less than approximately 2.5_⊙ the Roche radius is simply too small for the size of the star. We therefore conclude that the observations are inconsistent with a stripped star nature for the primary.

Fig. 13. Mean value of the Roche radius of the primary versus primary mass for the binary parameters relevant here. The area within the dashed lines (lower left) is prohibited due to the primary’s Roche radius shrinking to less than the size of the star. We plot the mean value of the Roche radius in this figure, however for all values of sin i greater than 0.1 the actual value is only a weak function of sin i.

4.2. The primary as a massive star

Adopting m₁ = M_ev = 8.6 M_⊙ for the primary we construct the diagnostic diagram in Fig. 14. Here we plot m₂ (the secondary mass), Roche₁ (the Roche radius of m₁) and the predicted f_el as a function of sin i. Also indicated, on the vertical axis, is our estimate of the radius of the primary, R₁ = 5.6 R_⊙, that fits comfortably inside its Roche radius, which is roughly twice the stellar radius, for all reasonable values of sin i. We use the analytic approximation of Morris (1985) for f_el and ignore the brightness of the secondary, that is, we assume that only the primary is bright. We also used the binary simulation package nightfall⁴ to model binary systems for specific main sequence secondary stars (see for example Table C.1), at appropriate values of sini. From the resulting light curves we also plot the maximum percentage variation we would expect to see in the data. Note that at high sini values the light curves should display eclipses, that are not seen in the data. The observed limits on f_el, labeled f_MACHO and f_Gaia for the MACHO and Gaia data respectively, are also indicated. The latter is particularly constraining, implying sin i <  0.4.

Fig. 14. Diagnostic diagram assuming a main sequence primary with a mass of 8.6 M⊙. The solid curves represent, as a function of sini, the secondary mass (labeled m2), the Roche radius of the primary (Roche1), and the maximum percentage amplitude of ellipsoidal light variations expected (fel) of the primary. Filled stars (⋆) illustrate the percentage amplitude of ellipsoidal light variations predicted by nightfall assuming a binary consisting of the proposed primary and a main sequence secondary of appropriate mass and radius. The increase at high sini indicates the presence of an eclipse. Dotted horizontal lines illustrate the constraints implied by the observations, the stellar radius (R1), the limits on ellipsoidal light variations (fMACHO and fGaia), and the maximum mass that is consistent with the absence of secondary stellar spectrum in the data. The vertical line (labeled Synch) indicates the upper limit on sin i implied by synchronous rotation and our measured upper limit on ve sin i of 15 km s−1.

The upper limit on sini, implied by synchronous rotation, is also indicated in Fig. 14, and moreover implies that sin i <  0.15. Hence the secondary should have a mass of at least m₂ ∼ 25 M_⊙, and should be a BH since we would easily see such a massive star. Normal stars in a binary of this period have circularization (τ_circ) and synchronization (τ_synch) timescales on the order of 10⁶ and 10⁴–10⁵ years respectively (Hurley et al. 2002), and for this specific case we estimate τ_synch = 0.07 Myr. We are led to the conclusion that the B-type star is likely in synchronous rotation around a stellar mass BH.

We add one additional constraint on the vertical axis, namely the mass ($m_2^{\text{max}}$) above which a narrow-lined main sequence secondary star would be visible in the spectrum. In Appendix C we discuss in detail how we make this estimate, finding that $m_2^{\mathrm{m}\mathrm{a}\mathrm{x}}=3.3M_{\odot }$. Ignoring for the moment the assumption of synchronous rotation, viable solutions for a slowly rotating secondary would imply a secondary mass range of 2.1–3.3 M_⊙. A fast rotating asynchronous secondary must be at least dimmer and less massive than the tertiary, but as is discussed in Sect. 4.3, the disk emission appears incompatible with the radius of such a star.

4.3. Hα and IR emission

As discussed in Sect. 3, it is unlikely that there is a classical Be star in the system as the Balmer emission appears to be associated with the inner binary, although it does not share the primary’s orbital motion. We therefore assume that the emission is centered on a more massive secondary. The morphology of the Balmer emission lines suggests a disk, and we can therefore use velocity measurements to estimate its size. The observational quantities from the VLT-FLAMES data from Fig. 10 are the peak and wing velocities of V_p = 155 km s⁻¹ and V_w = 500 km s⁻¹. The latter is of course a lower limit as it only measures the extent of the emitting disk visible in those data. While the SALT-HRS emission lines also indicate a possible double peaked structure, see Fig. 11, this is difficult to measure, hence we instead use the half-width half-maximum (HWHM) velocity as a characteristic property of the disk, with V_h = 27.5 km s⁻¹. On the assumption that these velocities reflect Keplerian motion in a disk we can infer some approximate dimensions for given values of i and m₂ by estimating the radii implied by these velocities. Roughly speaking V_p is related to the outer boundary radius, and V_w to the inner boundary radius (Huang 1972; Horne & Marsh 1986).

Our estimates of the radii assume that the measured velocities, corrected for sin i, match the Roche potential along the line joining the two masses (in the co-rotating frame, assuming synchronous rotation), that is, in the direction of the L1 Lagrange point. In reality, orbits about m₂ in the plane of rotation are elongated along this axis but comparison with numerical solutions (see for example Paczynski & Rudak 1980; Artymowicz & Lubow 1994) indicates this approach gives a reasonable estimate of the mean radius of stable orbits. This estimate improves with decreasing radius, though of course makes no allowance for orbit instabilities or tidal effects, that are more important at larger radii. Figure 15 illustrates the expected radii for the three characteristic velocities above, assuming a primary mass of 8.6 M_⊙, compared with the binary separation and the Roche radius of the secondary. We continue to refer to m₂ as the secondary even though we expect m₂ >  m₁ in the current scenario. We also show the locus of the secondary’s Roche radius, and the binary separation (a), together with curves representing the maximum disk radius permitted by tidal forces (Artymowicz & Lubow 1994), and maximum radius permitted by stability considerations (taken from Paczynski & Rudak 1980) of free-particle orbits. Similar results are found from more complex calculations for viscous disks by Artymowicz & Lubow (1994), though differences are found for nonzero eccentricity and increasing Reynolds number.

Fig. 15. Keplerian radii for characteristic Hα emission line velocities (dashed lines) as discussed in the text and labeled as follows; Vw is the maximum extent of the emission in the FSMS data, Vp is the velocity of that double peaked structure, and Vh is the HWHM of the SALT/HRS emission. The solid curves represent the binary separation (labeled a) and the Roche radius of the unseen secondary (Roche2). The dotted line is the limit imposed by tidal forces, while the dash-dotted line joining stars is the limiting average radius of stable orbits of the outer disk, from Paczynski & Rudak (1980).

Figure 15 suggests that the outer disk velocity V_p implies the disk can fit inside the Roche lobe provided sin i < 0.25. (The discontinuity in the V_p curve near this point is simply an artifact, due to the velocity exceeding that at the L1 Lagrange point. Values above this point should be ignored.) V_w on the other hand is well within the Roche radius for all sin i. The Hα emission profile from the VLT-FLAMES data therefore suggests the presence a disk at small values of sin i, roughly lying between 0.7 and 7 R_⊙, consistent with the results from the previous subsection.

Considering now V_h, from the SALT-HRS data, its value implies the gas lies at a very large distance from the binary, from a few hundred to a few thousand solar radii. Furthermore, the variability of the disk emission confirms that the proposed disk is unstable on timescales of years. We have no explanation for this other than to suggest the inner binary has somehow ejected its disk, possibly as a result of a close encounter with the tertiary.

The compactness of the disk, as implied by the VLT-FLAMES data, is consistent with the very faint IR excess observed in the SED, with an appropriate choice of uncertain disk parameters, see Fig. 16. For example a flat Be-star disk model (Carciofi & Bjorkman 2006) reproduces the observed IR excess for an inner disk temperature of ∼30 000 ± 10 000 K and inner/outer disk radii of 0.7/7.0 R_⊙. A simple α-disk accretion model (Shakura & Sunyaev 1973; Beall et al. 1984) requires an accretion rate of 10⁻⁷ M_⊙ yr⁻¹ in order to reproduce a similar IR excess, with a BH mass of 25 M_⊙. Truncating the outer radius at 7 R_⊙ reduces the IR excess significantly, hence the high accretion rate needed to match the observations. The latter would imply an X-ray luminosity several orders of magnitude in excess of the observed limit from XMM-Newton of L_x <  5.3  ×  10³³ erg s⁻¹, and is also far above the expected accretion rate of such a system. For example, a B-type star with the primary’s parameters is expected to have a mass-loss rate of ∼10⁻¹⁰ M_⊙ yr⁻¹ (Krtička 2014) giving an accretion rate onto the BH roughly an order of magnitude smaller than this. Decreasing the accretion rate in the model to this magnitude would require increasing the BH mass by some orders of magnitude to recover the IR excess. However, at such a low accretion rate the accretion may proceed via an advection dominated accretion flow (ADAF) as is the case for BH X-ray binaries at low accretion rates (Esin et al. 1997; Narayan & McClintock 2008). While producing ADAF models is beyond the scope of this paper, and in any case we only have upper limits on the X-ray flux, such models might succeed in matching the small IR excess, without exceeding the X-ray limiting flux. In summary, while the weak IR excess can be modeled with a disk size within the suggested limits, the nature of that disk is still obscure.

Fig. 16. Comparison of model SEDs (continuous lines) with extinction corrected photometric data (stars with labels, error bars are small and difficult to see on this scale). Models are normalized to the V-band and parameters are adjusted to fit the IR excess as noted in Sect. 4.3. Color coding as follows: black – composite stellar spectrum; magenta – flat Be disk and accretion disk model (both models are essentially identical on this scale); blue – for comparison this illustrates a flat Be disk model if it were circumstellar to the B-type star and truncated at its Roche radius.

Implicit in the above discussion is the assumption that the disk is in equilibrium. From the observations this is clearly not the case, though the strong Balmer emission in the VLT-FLAMES data is present for ∼60 days, or roughly 30 periods of the inner binary. However, considering circumbinary disks, binaries tend to clear the inner disk to a distance of ∼2a on dynamical timescales (Artymowicz & Lubow 1994; Rudak & Paczynski 1981), that is, within a few tens of periods. It is possible, therefore, that the disk is episodic and the emission observed at later times, in the SALT-HRS data, is a remnant circumbinary structure. We assume that the tertiary is essentially a spectator in this scenario.

4.4. Constraints from the tertiary

Since our knowledge of the tertiary star’s properties are only approximate we can only make some broad inferences. The V-band flux ratio yields an absolute magnitude of M_V = −2.5 m for the tertiary and, assuming it lies on an isochrone with the narrow-lined B-type star, implies a mass of approximately 6–7 M_⊙. While the longer period is very poorly sampled, the VLT-FLAMES data imply that a period of less than ∼120 d is inconsistent with the data. If we assume that spin and orbital angular momentum of the tertiary are aligned then its high v_e sin i argues for a significant tilt of its orbital plane with respect to the inner binary. Assuming an equatorial rotational velocity close to the escape velocity of the star provides lower limit of sin i > 0.5.

The stability of such a triple star system can be investigated following Toonen et al. (2020) and is illustrated in Fig. 17⁵. For these tests we adopted an inclination angle between the inner and outer binaries of 87° in order to show the limit of secular instability. Changing the inclination such that it differs from 90°, the angle that maximizes the secular limit extent, by only a few degrees more shrinks the secular limit to negligible size. We see that stability increases with increasing outer period and decreasing eccentricity and, from the above, that the limit of secular instability is only significant for inclinations very close to 90°.

Fig. 17. Visualization of stability criteria for the triple systems (Toonen et al. 2020, where further details can be found), the shaded area below the red line indicating unstable configurations. The locations of test models are illustrated, assuming m1 = 25, m2 = 8.6, m3 = 6 M⊙, and Pin = 2.92 d. The range of outer periods and eccentricities are indicated by the labels, for example the model p120ep1 assumes Pout = 120 d and an eccentricity of 0.1 for the outer system. All test models assume prograde orbits of 87° inclination of the orbital planes.

It is difficult to be more specific except to note that the velocity dispersion of the radial velocity measurements of the proposed tertiary are larger than those of the emission line velocities, and the possible range of inner binary systemic velocities (see Fig. 4 and Tables 4 and 5). Naively, one might interpret this as support for the BH hypothesis, though there are significant uncertainties for the tertiary results.

Further radial velocity monitoring would help clarify the nature of the system. For example, the primary has such narrow metal lines that radial velocity measurement is straightforward. Hence spectroscopic monitoring with a cadence covering the inner period, and an appropriate overarching cadence to sample potential outer orbital periods would enable two complementary methods. On the one hand, the radial velocities of the tertiary may be retrievable with sufficient accuracy to determine the mass ratio of the inner binary to the tertiary directly, thus enabling a determination of m₂. On the other hand, should that not be feasible, the systemic radial velocity data of the inner binary alone may permit determination of m₂ provided one can determine the mass function (f) for the system consisting of the inner binary and the tertiary, and hence solve the resulting quadratic equation to give

$$\begin{array}{c}\hfill m_1+m_2=\sqrt{m_3^3{\sin }^{3}i/f}-m_3\end{array}$$

subject of course to eccentricity and orbital plane considerations. An interesting corollary is that triple, or higher order, systems may provide an unambiguous method for determining the masses of suspected dark or subluminous members, depending on their periods.

5. Conclusions and caveats

We determine NGC 2004#115 to be a triple system comprising an inner short-period (2.92 d) system with a very low v_e sin i 8.6 M_⊙ B-type primary star and, assuming synchronous rotation, a BH of mass > 25 M_⊙. The tertiary is determined to be a broad-lined 6–7 M_⊙ B-type star with an outer period in excess of 120 days. No sign of a secondary spectrum is observed in the inner binary. The velocity structure of the early time doublepeaked Hα emission line is consistent with an accretion disk inside the Roche lobe of the secondary, and its centroid velocity is remarkably constant. It is therefore useful to consider the plausibility of forming an inner binary consisting of an 8.6 M_⊙ main sequence stars and a 25 M_⊙ BH.

While the current system is relatively stable (from Fig. 17) we have no knowledge of its prior state. Nevertheless, Toonen et al. (2022) recently presented an overview of the evolution of triple systems and showed that, ignoring the small percentage of disrupted systems, roughly three quarters of triples follow a path in which the primary star of the inner binary initiates mass transfer. The remaining evolve with no interaction except for a small percentage that follow paths in which either the secondary or tertiary initiates mass transfer, each comprising of roughly 1% of the total. We thus consider the evolution of the inner binary in isolation.

Assuming the BH mass equals the helium core mass at collapse implies a progenitor star of at least 50 _⊙ (Woosley 2019). This makes an initially short orbital period (≲20 d) rather unlikely, since most systems with such extreme initial mass ratios (< 0.2) cannot undergo stable mass transfer, and would therefore quickly merge (cf Fig. B1 in Langer et al. 2020a). Binaries that are initially wider may undergo stable mass transfer for extreme mass ratios if the mass transfer is inefficient; however, the post-mass transfer orbital periods would generally be larger than 10 d. Furthermore, the normal surface abundances of the inner B-type star argue against significant mass transfer from the primary.

Common envelope (CE) evolution may provide a more viable explanation. Stars of 50 M_⊙ or more are expected to go through an LBV phase (Smith et al. 2011), and lose their hydrogen-rich envelope even as single stars. Therefore it is plausible that an 8.6 M_⊙ companion in an orbit at several 100 R_⊙ separation could survive CE evolution, but be dragged into a short period orbit (Kruckow et al. 2016). Alternatively, the inner binary could result from a tight extreme mass ratio binary in which only the massive component evolves chemically homogeneously. According to Marchant et al. (2017), while this scenario requires slightly sub-LMC metallicity, it may become possible by enhanced mixing in the closest binaries (Hastings et al. 2020; Koenigsberger et al. 2021), which might also explain the shortest period Wolf-Rayet binaries in the LMC (Shenar et al. 2020b).

An important caveat is that the inclination hinges on the assumption that the inner binary is in synchronous rotation, the upper limit on sin i of 0.15 (9°) following from the period of the system and the radius of the primary B-type star. The period is well determined, in spite of the incomplete phase coverage of the VLT-FLAMES data, and leaves little room for manoeuvre. The radius is determined by the assumption that the star must have LMC composition and the system resides in the LMC. The latter condition is fully consistent with its negligible parallax, proper motion, systemic radial velocity and presence of LMC-like interstellar lines. A substantially smaller radius would imply significant helium enhancement (Lennon et al. 2021) accompanied by a very strong C/N depletion/enhancement. The latter is not supported by the data. Furthermore, as the cirularization time of the inner binary, ∼16 Myr, is much larger than its synchronization time, 0.07 Myr, this strongly supports the assumption of synchronicity. Indeed an 8.6 M_⊙ star reaches a radius of 5.6 R⊙ after about 20 Myr (H-exhaustion occurs after 32 Myr, Brott et al. 2011). This picture would be consistent with the BH-progenitor initiating interaction after about 4 Myr, the lifetime of a massive star.

Nevertheless, an inclination angle of 9° or lower has a probability of about 1%, hence making the ad hoc assumption that the inner binary is asynchronous, there are two options: If the star is rotating more quickly then i must be smaller than 9°, the BH mass more massive, and we are entering into an even less probable scenario. If the star is rotating more slowly, i is larger and the BH mass is smaller. Indeed if i is large enough we enter the range where the companion might even be a more mundane intermediate mass star. In this case, the primary would have a very low equatorial velocity and hence would need to lose most of its angular momentum. The most feasible way of achieving this is via Roche lobe overflow, leading to the formation of a stripped star. However, as already demonstrated, this solution is not viable. Furthermore, if the secondary would indeed be a main sequence star of similar mass to the component identified as the tertiary, its Balmer line radial velocity variations would be substantially larger than those already identified for the tertiary.

In conclusion, while our investigation into the nature of this enigmatic system leads to the surprising result that it may harbor a ∼25 M_⊙ black hole, further spectroscopic monitoring is required on cadences that will provide stronger constraints on the both inner binary and the proposed tertiary. Such data might also be amenable to spectral disentangling techniques.

During the refereeing stage of this paper we became aware that Saracino et al. (2022) reported the discovery of a B+BH binary in the LMC cluster NGC 1850, consisting of a 5 M_⊙ B-type star plus an 11 M_⊙ BH. However El-Badry & Burdge (2022) have argued that the primary star is a stripped star with a mass of only ∼1 M_⊙, possibly with a main sequence companion. Applying our methodology (cf. Fig. 14) to the B+BH system we estimate a stellar radius of 6.6 R_⊙ and a stellar luminosity of log L/L_⊙ = 3.2. The mass function yields a Roche radius of approximately 10 R_⊙, giving a fill factor of only 0.66. This system cannot reproduce the observed light curve as it would generate maximum (i.e., sin i = 1) ellipsoidal light variations with an amplitude of only ∼4–5%, significantly smaller than the observed light curve constraint of ∼10–12%, in broad agreement with El-Badry & Burdge (2022). Reducing the mass of the primary we find that its Roche radius matches the stellar radius for a mass of ∼1.8 M_⊙, the light curve then implying i  ∼  47° with a BH mass of ∼5.8 M_⊙. However if a secondary contributes significantly to the light then the radius and luminosity of the primary would be smaller, and the intrinsic ellipsoidal light variations would be larger, thus permitting an even lower mass primary, and secondary, as a possible solution. In this case the mass and luminosity of the primary would imply that it is a stripped star, as also found by El-Badry & Burdge (2022). As described in the present paper, a good quality high resolution spectrum may enable a detailed abundance analysis and a direct determination of the flux contribution and radius of the narrow-lined B-type component.

Acknowledgments

DJL acknowledges support from the Spanish Government Ministerio de Ciencia, Innovación y Universidades through grants PGC-2018-091 3741-B-C22 and from the Canarian Agency for Research, Innovation and Information Society (ACIISI), of the Canary Islands Government, and the European Regional Development Fund (ERDF), under grant with reference ProID2017010115. DJL thanks Ian Howarth for providing high resolution ATLAS spectra for the LMC from Howarth (2011) that were instrumental in looking for cool companions in the inner binary. The SALT data were acquired under the auspices of Director’s Discretionary Time proposal 2020-1-DDT-007 (PI: I. Monageng). We thank the referee, Tomer Shenar, for a careful reading of the manuscript that led to a significant improvement of the paper. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of the SIMBAD and VIZIER databases (operated at CDS, Strasbourg, France) and the table manipulation software TOPCAT (Taylor 2005).

Note added in proof: Using OGLE data El-Badry et al. (2022) have since shown that this system consists of three luminous stars, with the close inner binary being an asynchronous pair of stars with a mass ratio of approximately 2. They have also pointed out an error in our analysis of the Gaia data leading to an underestimate of the system’s photometric variability by a factor of 21.

References

Appendix A: Radial velocity measurements

A.1. VLT-FLAMES data

Radial velocities of the narrow lined B-type star were measured by cross-correlating each spectrum against an LMC model template of similar spectral type (an exact match is not required), convolved with the instrumental resolution (from Evans et al. 2006) and appropriate v_e sin i. Where possible we used the metal lines to measure the radial velocities, selecting segments of the spectrum in each of the observed FLAMES settings that contain sufficient strong metal lines to produce a strong correlation peak, avoiding the hydrogen and helium lines that may be broad, asymmetric or have a contribution from the disk. Note that He I lines were used for the HR02 observations due to the lack of suitable strong metal lines therein. The cross-correlation segment regions are listed in Table A.1, together with a summary of the metal species included in those segments. As noted in Sect. 2, each of the six regions was observed in two observing blocks (OB), OB1 and OB2, with three separate exposures in each OB, exposures 1–3, and 4–6. Radial velocities were measured for each individual exposure except for the HR03 observations as exposures 5 and 6 of OB2 had very low signal levels and were not useful. The blue spectrum, illustrated in Fig. A.1, is a typical sharp lined B2 spectrum, the B2e classification coming from the peculiar Hα emission line shown in Fig. 1.

Fig. A.1. VLT-FLAMES merged data for the spectral regions, top to bottom, HR02, HR03, HR04, HR05 and HR06. Approximate normalization of the spectra was accomplished with a low order polynomial, for line profile fitting a custom normalized was performed.

Radial velocities were measured by approximating the cross-correlation peak with a Gaussian function, a good approximation for the narrow metal lines, and determining its position using a least squared fit to its width, position and base level. The 1-σ uncertainty in the position is adopted as the uncertainty in the radial velocity. The extreme narrowness of the lines leads to typical uncertainties of less than 1 km s⁻¹ for each spectrum, though the systematic velocity precision from this instrument is estimated as ∼2 km s⁻¹ between grating settings (Evans et al. 2008). The results are given in Table A.2 and, depending on epoch, the short period enables one to discern small velocity shifts within the separate exposures of certain OBs.

As a check of these values we also determined the maximum of the peak position of the cross-correlation function as a measure of the radial velocity. Both methods were found to give very similar results, with agreement within the errors, except for the HR02 data where small shifts of up to 5 km s⁻¹ are observed between these two methods that we attribute to the use of He I lines that are slightly asymmetric and skew the Gaussian fit away from the peak position.

A.2. SALT-HRS data

The HRS instrument on SALT is a dual blue-arm/red-arm echelle spectrograph that was used to cover the spectral regions 3730–5570 and 5420–8790Å at a spectral resolution of R∼ 14 000. The S/N of the observations for the blue arm data ranged from 20–40 per pixel, though it dropped very rapidly below 4500Å such that the lines in this region were not useful. The S/N for the red arm data was more uniform and ranged from 20–30 per pixel, although the night sky lines were very strong. This latter issue was a concern for the Hα region and it was decided not to subtract the sky lines from the data to avoid unknowingly distorting the profile. Due to the low S/N of the SALT data radial velocities were determined using a Gaussian fitting approach to lines belonging to the Si III triplet at 4553/4568/4575 Å and the C II doublet at 6578/6383 Å. The low S/N in the blue region precluded the use of metal lines below 4500Å. Results are also listed in Table A.2.

Appendix B: Atmospheric analysis of the B-type spectrum

Typically, T_eff and log g are determined by satisfying one or more of the ionization balances Si II/Si III, Si III/Si IV or He I/He II, while simultaneously fitting Balmer lines such as the Hγ and Hδ profiles. In the case of NGC 2004#115, Si III is the only ionization species that is present with significant strength, both Si II and Si IV are marginally present, while He II is absent. Nevertheless the weakness of the Si II and Si IV lines provide strong constraints on the effective temperature, but the Balmer lines are compromised due to the presence of another component preventing their use as the main log g constraint.

Therefore, for a given value of log g we determined the effective temperature and element abundances by comparing observed equivalent widths of the metal lines, Table B.1, with precomputed grids. Table B.1 includes the C II doublet at 6580Åthat was not included in previous FSMS analyses, but has been found by Barnett & McKeith (1988) to provide reliable carbon abundance estimates. Excluding this doublet would change the estimates by less than 0.05 dex. The veiling effect of second light is approximated by a stellar continuum with T_eff = 20 000 K although changing this temperature by ±5 000 K has little impact on the final results. This continuum is scaled to the required V-band flux ratio and used to rescale the predicted equivalent widths to diluted values for comparison with observation. The overarching procedure then consists of running a sequence of calculations at a given log g to find the T_eff and flux ratio that reproduces the normal LMC silicon abundance, the results of which are presented in Fig. 5. Further details on the method can be found in Lennon et al. (2021).

The microturbulent velocity (v_t) deserves special mention. It was estimated from the relative strength of the components of the Si III triplet near 4560Å, a positive abundance gradient with equivalent width indicating too low a value, a negative gradient, too high a value. With no veiling, a negative slope is found for all values of the v_t, including v_t = 0 km s⁻¹. This anomaly has been identified previously by Hunter et al. (2007) for a small number of FSMS targets, with one explanation being veiling by an unidentified companion. For larger veilings (∼50%), the observed equivalent widths are consistent with v_t ≃ 0 within their estimated errors. Hence we conclude that the microturbulent velocity is small with our best estimate being 0 km s⁻¹, though with low significance as values up to 2–3 km s⁻¹ are difficult to exclude. The main impact of an increased v_t is to increase the veiling (i.e., decrease the contribution of the B-type star to the total flux), for example with v_t = 2 km s⁻¹ the B-type star would be fainter by about 6% of the total flux.

The projected rotational velocity, v_e sin i, was estimated from the Si III 4560Å triplet which was observed with both the HR05 and HR06 settings, providing six independent estimates. Theoretical spectra for the best fitting atmospheric parameters listed in Table 2 were convolved with the instrumental (assumed to have a Gaussian profile and with FWHMs taken from Evans et al. 2006) and rotational profiles, and then re-binned to match the observations. In general, the observed profiles are best fitted with v_e sin i∼10 km s⁻¹, with the estimates being little affected by changes in the adopted effective temperature or gravity. Increasing the microturbulent velocity to 2 km s⁻¹ would have reduced the estimates by typically 2 km s⁻¹. Similarly changing the FWHM of the instrumental profile by 10% would typically change the v_e sin i estimates by 1 km s⁻¹. In summary our best estimate of the projected rotational velocity is 10 km s⁻¹. Allowing for possible systematic uncertainties leads to a robust upper limit, v_e sin i ≤ 15 km s⁻¹.

Appendix C: Limits on a secondary spectrum

Here we examine whether the secondary, if it is a main sequence star, would be visible in our data. We make the assumption that this secondary is also narrow lined, by virtue of the expected short synchronization times. Since a narrow-lined secondary is not observed in our data this implies that the radius of the secondary is smaller than the primary, and its v_e sin i should be smaller. Hence, its spectrum would not be resolved with the present data and we therefore assume a line width determined by instrumental resolution of FWHM∼15 km s⁻¹. A broad lined, or quickly rotating, star would clearly be much more difficult to detect. We adopt notional effective temperatures and luminosities for a range of near main-sequence stars, Table C.1, and assign them spectra from the TLUSTY NLTE grid discussed in Appendix B. Expected masses are also listed in this table for two different age assumptions, where we have used the zero rotation isochrones of Mowlavi et al. (2012), as the tracks of Brott et al. (2011) do not extend to masses below 5 M_⊙. Note that evolution has little impact for the fainter stars, the same being true for rotation, and indeed the exact choice of stellar evolution models.

As more massive and hotter stars would be easily identified due to their contribution to the hydrogen Balmer lines and He I lines (Fig. C.1), we focus on features that provide a strong contrast for the cooler stars. The Mg II doublet at 4481 Å is well suited to this purpose. It is a strong feature that has a maximum strength on the main sequence at roughly A0 spectral type, hence strengthening as one moves to cooler, less massive stars. There are other candidates, such as the Si II doublet at 6347,6371 Å and the Ca II doublet at 3934 Å (the 3968 Å component is blended with Hϵ). However, Mg II is preferred as the signal-to-noise in this region is superior, it lies in a relatively clean part of the spectrum, and moreover the relative radial velocities of OB1 and OB2 are significantly different from zero; approximately -28 and +80 km s⁻¹ respectively relative to the systemic velocity (values are adopted from the orbital solution). The negative shift of the primary is particularly useful, as the proposed tertiary is also blue shifted at these epochs, hence the secondary should be red-shifted at that time.

Fig. C.1. Comparison of observed data for HR05 OB1 and OB2 (black lines) with simulations of the sum of primary, secondary and tertiary spectra (blue lines) weighted according to their absolute visual magnitude, MV, from Table C.1. The main features are the He I 4471 Å and Mg II 4481 Å lines. The former shows no evidence of the secondary implying that the spectral type must be B3 V or later. The constraint derived from the Mg II is discussed in the text.

Synthetic secondary spectra were computed assuming a radial velocity appropriate for the assumed mass, weighted by their respective magnitudes, corrected to distance of the LMC, and then diluted with a featureless continuum representing the tertiary plus primary such that the observed dereddened magnitude of the system is conserved. The equivalent width of the resulting Mg II line of the secondary was then measured and, since no such component has been identified, it was compared with the detection limit (limiting equivalent width) of the data. Since the secondary component has different positions in OB1 and OB2 we shifted both OBs to the rest frame of the potential secondary, then merged the spectra and measured the limiting equivalent width in the resulting spectrum. The results are summarized in Table C.1 and can be compared with the 1-σ detection limit in this region of the spectrum that we estimate as 2.6 mÅ (based on 1000 equivalent width measurements at random continuum positions). We can therefore exclude, at a 2-σ confidence level, the presence of a secondary main sequence star with an absolute visual magnitude above ∼+0.6. If the star would be a zero age main sequence star its mass would be ∼3.3 M_⊙, a more evolved star would be correspondingly less massive.

Allowing higher values of v_e sin i for the secondary obviously makes it harder to detect. Values of 100 and 200 km s⁻¹ imply 1-σ detection limits of 6.6 and 8.8 mÅ respectively for this line. On the other hand the brighter stars, B3 V and above, are more easily detected using their stronger H I and He I lines for these moderate v_e sin i values.

Appendix D: Spectral energy distribution

Sources of optical photometry listed in Table D.1 include the MCPS catalog of Zaritsky et al. (2004), Evans et al. (2006, Evans) and Gaia (Gaia Collaboration et al. 2018). This source lies within the boundaries of the MACHO survey of the LMC (Alcock et al. 2000) and was monitored for a period of almost 7 years. Mean magnitudes in the Kron-Cousins system of V = 15.50 ± 0.04 and R = 15.64 ± 0.04 are derived from the instrumental magnitudes using the transformations and zero-points described in Alcock et al. (1999). Near-IR data are from the IRSF Magellanic Clouds catalog (Kato et al. 2007), which we prefer to 2MASS as the latter only contains upper limits for the H and K_s magnitudes, and the uncertainty of J is 0.1^m. The mid-IR magnitudes are from the Spitzer SAGE survey (Meixner et al. 2006), and are taken from the catalog of Bonanos et al. (2009). The source was not detected longwards of the 4.5μm band.

Stroud (2011) used the MACHO data to investigate the field of NGC 2004 in a search for variable stars, finding that the light curves of NGC 2004#115, in the MACHO V and R bands, are remarkably constant. A Lomb-Scargle test on both MACHO V and R bands does not find signs of the spectroscopic period. Two peaks above the 0.1% false-alarm-probability level are found at 3.2 and 2.2 days in the case of the V and R bands respectively. However, the discrepancy between these two periods and the shape of the periodograms suggest that these are spurious peaks, Fig. D.2. Folding the photometry data to the spectroscopic period (Fig. 12) we estimate 1-σ limiting peak-to-peak amplitudes of 0.002 mag, in the case of the V band, and 0.006 mag for the R band. These constrain the amplitude of possible ellipsoidal light variations to less than 0.005^m in V and 0.007^m in R in amplitude, or around 0.5–0.6%. The lack of variability of the system is also reflected in the Gaia EDR3 data. While the project has yet to release all the epoch photometry we note that its G-band flux is given as 1.2e4 e⁻ s⁻¹, with a 1-σ error of 10.73 e⁻ s⁻¹ based on 459 observations. In fact, according to Gaia data, NGC 2004#115 is one of the least variable B-type stars observed in NGC 2004 by the FSMS, as illustrated in Fig. D.1. While Hα displays significant profile variability, this has negligible impact on broad-band magnitudes.

Fig. D.1. GaiaG magnitude is plotted against its % flux error, as an indication of variability for the B-type stars observed in NGC 2004 by Evans et al. (2006). Star #115 (inverted triangle, lower left) is among the least variable stars in this plot. Candidate classical Be stars are circled (from Dunstall et al. 2011).

The field of this system was also the target of two pointed XMM-Newton observations of ∼5000 seconds each, on 2007-01-28 and 2011-12-20, although there are no obvious X-ray sources in either data-set near its position. We used a 15″radius aperture to obtain 2-σ upper limits of 0.0066 and 0.0034 counts/sec respectively in the 0.2–12 keV band. However, combining the images to get a total exposure time of 10.2 ks results in an upper limit of 0.0055 counts/s, as the higher background in the first observation dominates the combined image. Converting the lowest of these limits into flux using using a power-law of slope of two, and hydrogen column density of N_H = 2.7 × 10²¹ cm⁻², yields limits for the absorbed flux of 1.13 × 10⁻¹⁴ ergs cm⁻² s⁻¹, and for the unabsorbed flux, 1.78 × 10⁻¹⁴ ergs cm⁻² s⁻¹. This implies an X-ray flux from the system of less than 5.33 × 10³³ ergs s⁻¹ in the 0.2–12 keV band.

Given the possible composite nature of the SED we adopted the mean extinction for the cluster of E(B − V) = 0.08 (Caloi & Cassatella 1995; Cassatella et al. 1987) as appropriate for NGC 2004#115, and corrected the photometry for extinction using a combined average MW and LMC extinction curve (Gordon et al. 2003). Using the filter function taken from Bessell (1990), filter weighted model V-band fluxes are normalized to the observed V-band flux to yield the ratio of stellar radius to distance (R/d). From Table D.1 it is apparent there are some small zero point offsets in the various V magnitudes (more so in the B band), hence we simply adopt the MACHO result in determining the radius in sections 3 and 4. The SED, illustrated in Fig. 16, exhibits only a marginal IR excess at 4.5μm when compared with a model having no disk contribution.

Fig. D.2. Zoom of Lomb-Scargle periodograms of the MACHO V and R band photometry below 10 days illustrating the lack of a significant peak in either band.

All Tables

All Figures

Fig. 1. Composite nature of the Hα profile that consists of a broad double-peaked emission feature, plus a narrow absorption component from the B-type star, the narrow lines of the C II 6578/6582 Å doublet are further to the red. The two epochs from the VLT-FLAMES observations shown, OB1 (red) and OB2 (blue), illustrate significant velocity variability of the narrow components but a relatively stable broad Hα emission.

In the text

Fig. 2. Left panel: comparison of Hγ profiles from two VLT-FLAMES epochs in the rest frame of the narrow-lined B-type star. At these epochs the star is close to each quadrature phase. An arrow indicates the extra absorption that we attribute to the presence of a third star, not present in the blue spectrum. If this component was due to a faint companion to the B-type star then it should be visible on the red wing of the blue spectrum. Instead its velocity at that epoch is almost identical to that of the B-type star, see Sect. 3.2, hence no absorption is apparent at that epoch. Right panel: as for the left panel except here we plot two VLT-FLAMES epochs of the He I 4471 Å line. Similar considerations apply.

In the text

Fig. 3. Radial velocity (RV) orbital solution for NGC 2004#115 based on velocity estimates from the metal lines in its spectrum. Left-most panel: illustrates the fit using the VLT-FLAMES data only. Central panel: how the fit, now determined using both VLT-FLAMES (black points) and SALT-HRS (red points) data, deteriorates when we include both datasets. Rightmost panel: fit to the VLT-FLAMES data taken within the final week of that observing campaign, and fit parameters for this solution are listed in Table 1. The scales differ from panel to panel, in particular note the change in the scale of the O−C plots.

In the text

Fig. 4. Central panel: comparison of the solution derived from the late epoch FSMS data only (solid line) with the observed radial velocities of the primary, filled circles plus (small) error bars. Left panel: comparison of that solution with the early epoch FSMS data. Right panel; comparison with the SALT/HRS data. Horizontal lines indicate the systemic velocity of this solution. The dotted lines in the right panel are the solutions obtained by fitting those radial velocities but fixing both the period, eccentricity and semi-amplitude velocity to the values derived from the late FSMS data. Open star symbols are the velocities of the emission features (Table 5), and open circles are velocities of the broad-lined component (Table 4) that we identify as a possible tertiary.

In the text

	Fig. 5. Locus of (log g, Teff) parameters, labeled with the fractional B-star contributions, that reproduce the normal LMC silicon abundance of 7.2 dex. A helium abundance of N[He/H] = 0.1 and a microturbulence of vt = 0 km s−1 have been assumed.
In the text

Fig. 6. Subtracted spectra of the proposed tertiary for the HR02, HR04 and HR05 regions. We have not included HR03 due to its low signal-to-noise (discussed in the text), or HR06 as it is relatively featureless. Close examination of this spectrum may show reveal the presence of over-subtracted metal lines, but these are artifacts due to the use of a specific model spectrum from our grid with mean surface abundances a little different from the final abundance determination of the primary. Spectra were approximately normalized using a low order polynomial fitted to continuum regions.

In the text

Fig. 7. Illustration of the residual spectra obtained by subtracting theoretical spectra scaled by the fractional contribution from the HR06 OB2 data in the vicinity of He I 4471 Å line. The value of log g increases in steps if 0.25 dex from log g = 3.00 dex for the lower curve to 4.25 dex for the uppermost curve, spectra are offset by 0.2 normalized flux units. Teff values follow the relationship shown in Fig. 5. The emission feature at the position of the He I line core (dashed line), caused by over-subtraction of the Doppler core in the model, becomes weak or absent for log g > 3.75.

In the text

	Fig. 8. Residual (i.e., primary subtracted) Hβ profiles from the SALT/HRS data, observations 1 through 4 are from bottom to top, shifted by 0.5 units for clarity. The vertical lines indicate the mean positions of the absorption components (dotted line) and emission components (dashed line), from Tables 4 and 5.
In the text

Fig. 9. Hα profile morphology in the VLT-FLAMES data (top two profiles) and SALT-HRS data (four lower profiles). Spectra have been off-set by 0.5 normalized flux units for clarity. Sharp emission features in the SALT-HRS data are night sky lines (⊗). Hα and the sky line just long-ward of 6590 Å are likely blended with weak nebular H I and [N II] respectively, as discussed in the text. The position of the central dip of the Hα Doppler core of the narrow-lined B-type star is indicated, as are the positions of its two C II lines (↑).

In the text

Fig. 10. Left panel: VLT-FLAMES OB1 (blue) and OB2 (magenta) Hα profiles resulting from subtraction of the underlying spectrum of the narrow-lined B-type star. Horizontal lines, same color coding, represent the rest frame velocity of the primary at these epochs. Note that a small residual velocity shift is still apparent (cf Fig. 1). Right panel: the same profiles plotted in their rest frame velocities; OB1 and OB2 spectra in blue and magenta respectively, the merged spectrum is plotted in black. Also illustrated are the velocities of the double peaks (±155 km s−1), the full-width half-maximim (FWHM) of the emission (550 km s−1), and the extent to which emission is present in the line wings (±500 km s−1). These spectra were taken ∼34 days apart and some variability is apparent within the double-peaked structure of the profile.

In the text

Fig. 11. SALT-HRS line profile morphology in velocity space of Hα (upper four curves) and Hβ (lower four curves), color coded by observation; black – obs1, blue – obs2, magenta – obs3 and orange – obs4. The primary’s underlying profiles have been subtracted from these data. While the Hα profiles exhibit some variability of intensity their velocities are quite constant, and are well modeled with two-component Gaussian fit with velocities at 314 km s−1 (for the ambient nebular emission) and ∼325 km s−1 (for the source emission). The latter has a FWHM of ∼60 km s−1.

In the text

	Fig. 12. MACHO V and R photometry spanning MJD 49001.6091 to 51542.5763, folded to the 2.92 d period of the system, illustrating the lack of variability in both bands. The magnitudes are instrumental; zero points and transformations to Kron-Cousins V and R magnitudes can be found in Alcock et al. (1999).
In the text

Fig. 13. Mean value of the Roche radius of the primary versus primary mass for the binary parameters relevant here. The area within the dashed lines (lower left) is prohibited due to the primary’s Roche radius shrinking to less than the size of the star. We plot the mean value of the Roche radius in this figure, however for all values of sin i greater than 0.1 the actual value is only a weak function of sin i.

In the text

Fig. 14. Diagnostic diagram assuming a main sequence primary with a mass of 8.6 M⊙. The solid curves represent, as a function of sini, the secondary mass (labeled m2), the Roche radius of the primary (Roche1), and the maximum percentage amplitude of ellipsoidal light variations expected (fel) of the primary. Filled stars (⋆) illustrate the percentage amplitude of ellipsoidal light variations predicted by nightfall assuming a binary consisting of the proposed primary and a main sequence secondary of appropriate mass and radius. The increase at high sini indicates the presence of an eclipse. Dotted horizontal lines illustrate the constraints implied by the observations, the stellar radius (R1), the limits on ellipsoidal light variations (fMACHO and fGaia), and the maximum mass that is consistent with the absence of secondary stellar spectrum in the data. The vertical line (labeled Synch) indicates the upper limit on sin i implied by synchronous rotation and our measured upper limit on ve sin i of 15 km s−1.

In the text

Fig. 15. Keplerian radii for characteristic Hα emission line velocities (dashed lines) as discussed in the text and labeled as follows; Vw is the maximum extent of the emission in the FSMS data, Vp is the velocity of that double peaked structure, and Vh is the HWHM of the SALT/HRS emission. The solid curves represent the binary separation (labeled a) and the Roche radius of the unseen secondary (Roche2). The dotted line is the limit imposed by tidal forces, while the dash-dotted line joining stars is the limiting average radius of stable orbits of the outer disk, from Paczynski & Rudak (1980).

In the text

Fig. 16. Comparison of model SEDs (continuous lines) with extinction corrected photometric data (stars with labels, error bars are small and difficult to see on this scale). Models are normalized to the V-band and parameters are adjusted to fit the IR excess as noted in Sect. 4.3. Color coding as follows: black – composite stellar spectrum; magenta – flat Be disk and accretion disk model (both models are essentially identical on this scale); blue – for comparison this illustrates a flat Be disk model if it were circumstellar to the B-type star and truncated at its Roche radius.

In the text

Fig. 17. Visualization of stability criteria for the triple systems (Toonen et al. 2020, where further details can be found), the shaded area below the red line indicating unstable configurations. The locations of test models are illustrated, assuming m1 = 25, m2 = 8.6, m3 = 6 M⊙, and Pin = 2.92 d. The range of outer periods and eccentricities are indicated by the labels, for example the model p120ep1 assumes Pout = 120 d and an eccentricity of 0.1 for the outer system. All test models assume prograde orbits of 87° inclination of the orbital planes.

In the text

	Fig. A.1. VLT-FLAMES merged data for the spectral regions, top to bottom, HR02, HR03, HR04, HR05 and HR06. Approximate normalization of the spectra was accomplished with a low order polynomial, for line profile fitting a custom normalized was performed.
In the text

Fig. C.1. Comparison of observed data for HR05 OB1 and OB2 (black lines) with simulations of the sum of primary, secondary and tertiary spectra (blue lines) weighted according to their absolute visual magnitude, MV, from Table C.1. The main features are the He I 4471 Å and Mg II 4481 Å lines. The former shows no evidence of the secondary implying that the spectral type must be B3 V or later. The constraint derived from the Mg II is discussed in the text.

In the text

	Fig. D.1. GaiaG magnitude is plotted against its % flux error, as an indication of variability for the B-type stars observed in NGC 2004 by Evans et al. (2006). Star #115 (inverted triangle, lower left) is among the least variable stars in this plot. Candidate classical Be stars are circled (from Dunstall et al. 2011).
In the text

	Fig. D.2. Zoom of Lomb-Scargle periodograms of the MACHO V and R band photometry below 10 days illustrating the lack of a significant peak in either band.
In the text