Discovery and Characterization of a Rare Magnetic Hybrid β Cephei Slowly Pulsating B-type Star in an Eclipsing Binary in the Young Open Cluster NGC 6193

In order to obtain initial estimates of the component effective temperatures (T_eff ) and radii (R), we performed a multicomponent fit to the combined-light, broadband spectral energy distribution (SED) of the HD 149834 system, starting with the T_eff and R estimates for the primary star from the spectroscopic study of Huang & Gies (2006a). The resulting T_eff and R were iteratively updated based on the joint light-curve and radial-velocity model (Section 4.2) until a final satisfactory SED fit was produced.

As shown in Figure 7 (black curve), the SED from 0.3 μm to 10 μm can be very well fit by a single component with T_eff = 21467 ± 246 K, log g = 3.95 ± 0.15, and [Fe/H] ≈ 0 as reported from a previous spectroscopic analysis in the literature (Huang & Gies 2006b), with a best-fit extinction of A_V = 1.4 ± 0.1 (consistent with other determinations of the reddening of this cluster; e.g., Rangwal et al. 2017).

Zoom In Zoom Out Reset image size

Interestingly, there is an apparent excess in the UV. In principle this could be due to activity on the low-mass companion star; young low-mass stars have ubiquitous UV excess even after they stop accreting. Alternatively, the UV excess could suggest an additional hot source in the system. However, adding a second component cannot reproduce that excess if the same A_V is applied. For example, adopting the secondary star's radius determined from the light-curve analysis (Section 4.2), we can reproduce the UV excess with a secondary T_eff = 25,000 K (the dark blue curve in Figure 7); however, this requires zero extinction. Reddening it by the same A_V as above yields the light blue curve, which predicts a flux ratio in the TESS band of 6.9%, much too large for the ∼2.5% flux ratio that is required by the light-curve analysis (Section 4.2).

Therefore the UV excess is not easily explained, and we also note evidence for an excess at 20 μm. We are able to reproduce the ∼2.5% flux ratio in the TESS band with the same A_V applied if we instead adopt T_eff = 13,000 K (red curve in Figure 7). Such a faint companion would require extraordinarily strong intrinsic UV emission to produce the observed UV excess, some five orders of magnitude over photospheric (see Figure 7). Thus the origin of the UV emission, if real, remains unclear (but see Section 4.4 regarding observed X-ray emission in the system).

Most importantly, the SED analysis provides robust initial constraints on the properties of the stars that we can utilize in the full solution below (Section 4.2). Adopting the parallax estimated in Section 2 and correcting for the 0.029 mas systematic offset of the Gaia DR2 parallaxes reported by the Gaia mission, we obtain via the bolometric flux from the SED fit an estimate for the primary star radius: R₁ = 4.96 ± 0.13 R_⊙. This radius, together with the spectroscopic $\mathrm{log}g$, then gives an estimate for the primary mass of M₁ = 8.1 ± 1.3 M_⊙.

We used the detrended TESS photometry described in Section 3, along with its internal uncertainties, to perform a light-curve analysis of HD 149834 using the eb code of Irwin et al. (2011), which is based on the Nelson-Davis-Etzel binary model (Etzel 1981; Popper & Etzel 1981) underlying the popular EBOP program by those authors. This model is adequate for well-detached systems such as HD 149834, having nearly spherical stars. The main adjustable parameters we considered are the orbit period (P), a reference epoch of primary eclipse (T₀, which in this code is strictly the time of inferior conjunction), the central surface brightness ratio in the TESS bandpass (J ≡ J₂/J₁), the brightness level at the first quadrature (m₀), the sum of the relative radii normalized by the semimajor axis (r₁ + r₂), the cosine of the inclination angle ($\cos i$), and the eccentricity parameters $e\cos \omega$ and $e\sin \omega$, with e being the eccentricity and ω the longitude of periastron. The radius ratio (k ≡ r₂/r₁), normally also an adjusted parameter in light-curve solutions, was instead derived at each iteration of our analysis based on other information, as we explain below. In order to account for the flux contamination in the TESS aperture described in Section 3, we also solved for the third light parameter ℓ₃, defined such that ℓ₁ + ℓ₂ + ℓ₃ = 1, in which ℓ₁ and ℓ₂ for this normalization are taken to be the light at first quadrature.

A linear limb-darkening law was adopted for this work, with coefficients for solar metallicity taken from the tabulation by Claret (2017) for the TESS band, interpolated to the temperatures and surface gravities of the stars based on our final analysis. The same source was used for the gravity darkening coefficients. Due to the large difference in temperature between the components indicated earlier, there is a strong reflection effect in the light curve caused by the energy from the primary being re-radiated toward the observer by the cooler secondary star, with a strength controlled by the reflection coefficient of the secondary, η₂ (Etzel 1981). We therefore solved for this additional parameter, while leaving the corresponding coefficient for the primary set to zero, as its impact is negligible. We accounted for the finite time of integration of the TESS photometry by oversampling the model light curve and then integrating over the 30 minute duration of each cadence prior to the comparison with the observations (see Gilliland et al. 2010; Kipping 2010).

In addition to the photometry, we incorporated the radial velocities of Arnal et al. (1988) directly into the analysis together with their reported uncertainties, solving for the velocity semi-amplitude of the primary star (K₁) and the center-of-mass velocity of the system (γ). The use of these older data extends the time-base significantly to nearly 35 years, thereby improving the period determination.

The solution was carried out within a Markov Chain Monte Carlo (MCMC) framework using the emcee code of Foreman-Mackey et al. (2013), with 100 walkers of length 15,000 each, after discarding the burn-in. We adopted uniform priors over suitable ranges for most parameters mentioned above, and a log-uniform prior for η₂. The prior for ℓ₃ was assumed to be Gaussian, with a mean value of 4% as derived in Section 3, and a standard deviation of 1%. Convergence of the chains was checked visually, requiring also a Gelman–Rubin statistic of 1.05 or smaller for each parameter (Gelman & Rubin 1992). The relative weighting between the photometry and the radial velocities was handled by including additional adjustable parameters f_phot and f_RV to rescale the observational errors. These multiplicative scale factors were solved for self-consistently and simultaneously with the other orbital quantities (see Gregory 2005), adopting log-uniform priors for both.

As we only have radial velocities for the primary star, an estimate of the mass M₁ of that component is required before we can derive the absolute mass and radius of the secondary. The spectroscopic study of Huang & Gies (2006b) reported an estimate of the surface gravity for the primary of log g₁ = 3.955 ± 0.023, although they considered 0.15 dex to be a more realistic estimate of the uncertainty. We adopt this more conservative error here. The log g value together with the radius estimate R₁ for that star from our SED fit above yields the necessary information to infer M₁. In order to propagate the uncertainties more accurately to the radius and mass of the secondary, we chose to add log g₁ and R₁ as adjustable parameters in our solution, constrained by Gaussian priors given by the measured values and corresponding errors of those two quantities. At each iteration we then solved for the mass and radius of the secondary, which immediately provide the radius ratio k ≡ r₂/r₁ = R₂/R₁ needed to model the light curve.

The results of our analysis for HD 149834 are presented in Table 1, in which we report the mode of the posterior distribution for each parameter. Posterior distributions of the derived quantities listed in the bottom section of the table were constructed directly from the MCMCs of the adjustable parameters involved. Our adopted light-curve model along with the TESS observations can be seen in Figure 8, with residuals shown at the bottom of each panel.

Zoom In Zoom Out Reset image size

Table 1. Results from Our MCMC Analysis of HD 149834

Parameter	Value	Prior
P (days)	${4.595701}_{-0.000038}^{+0.000035}$	[4, 5]
T0 (HJD−2,400,000)	${58641.8547}_{-0.0011}^{+0.0010}$	[58641, 58642]
J	${0.21}_{-0.11}^{+0.23}$	[0.0, 1.0]
r1 + r2	${0.249}_{-0.016}^{+0.014}$	[0.01, 0.60]
$\cos i$	${0.224}_{-0.017}^{+0.015}$	[0, 1]
$e\cos \omega$	$-{0.0009}_{-0.0017}^{+0.0014}$	[−1, 1]
$e\sin \omega$	$-{0.037}_{-0.014}^{+0.012}$	[−1, 1]
m0 (mag)	${7.12210}_{-0.00020}^{+0.00014}$	[7.0, 7.5]
η2	${0.00777}_{-0.00011}^{+0.00012}$	[−15, −1]
ℓ3	${0.0399}_{-0.0095}^{+0.0103}$	G(0.04, 0.01)
γ (km s−1)	$-{34.7}_{-1.8}^{+1.7}$	[−50, −20]
K1 (km s−1)	${22.5}_{-2.2}^{+1.9}$	[10, 50]
fphot	${8.54}_{-0.15}^{+0.19}$	[−5, 5]
fRV	${2.17}_{-0.26}^{+0.57}$	[−5, 5]
R1 (R☉)	${4.96}_{-0.14}^{+0.12}$	G(4.96, 0.13)
log g1 (cgs)	${4.040}_{-0.029}^{+0.034}$	G(3.955, 0.15)
Derived quantities
r1	${0.1940}_{-0.0054}^{+0.0035}$	⋯
r2	${0.0468}_{-0.0065}^{+0.0231}$	⋯
k ≡ r2/r1	${0.235}_{-0.026}^{+0.129}$	⋯
i (degrees)	${77.07}_{-0.86}^{+1.02}$	⋯
e	${0.037}_{-0.014}^{+0.012}$	⋯
ω (degrees)	${91.0}_{-1.8}^{+3.7}$	⋯
K2 (km s−1)	${251.6}_{-7.4}^{+11.0}$	⋯
M1 (M☉)	${9.71}_{-0.81}^{+1.22}$	⋯
M2 (M☉)	${0.858}_{-0.096}^{+0.123}$	⋯
q ≡ M2/M1	${0.0874}_{-0.0079}^{+0.0093}$	⋯
R2 (R☉)	${1.16}_{-0.13}^{+0.64}$	⋯
log g2 (cgs)	${4.11}_{-0.24}^{+0.23}$	⋯
a (R☉)	${25.57}_{-0.76}^{+0.97}$	⋯

Note. The values listed correspond to the mode of the respective posterior distributions, and the uncertainties represent the 68.3% credible intervals. All priors are uniform over the specified ranges, except those for η₂, f_phot, and f_RV, which are log-uniform, and the priors for R₁, log g₁, and ℓ₃, which are Gaussian with mean and standard deviations as indicated by the notation G(mean, σ).

Download table as:  ASCIITypeset image

Due to the reflection effect, the light ratio between the components of HD 149834 changes by a factor of about 1.9 as a function of orbital phase. This is illustrated in Figure 9. Our model radial-velocity curve with the observations of Arnal et al. (1988) is shown in Figure 10, indicating a center-of-mass velocity consistent with membership in NGC 6193, as mentioned earlier.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 11 depicts the two components of the HD 149834 system in the mass–radius and T_eff–radius planes, compared to the PAdova and TRieste Stellar Evolution Code (PARSEC) model isochrones (Bressan et al. 2012). In the mass–radius plane, the secondary appears roughly as expected for the nominal cluster age of 5 Myr, whereas the primary appears too large for that age; it appears slightly evolved with an inferred age of ∼20 Myr.

Zoom In Zoom Out Reset image size

In the T_eff–radius plane, both components can be reasonably well fit by the 20 Myr isochrone (note the very large uncertainty on the secondary T_eff due to the large uncertainty on J from the light-curve analysis). However, this is misleading because the secondary's measured mass of ∼0.86 M_⊙ should place it at a much cooler T_eff.

It is possible to explain both stars in the mass–radius diagram at about 5 Myr with a super-solar metallicity; however, the required metallicity is quite high at +0.7 (red curve). Unfortunately, there is no reliable literature estimate of the metallicity for HD 149834. Nonetheless, given all of the available evidence, we prefer this solution, as both stars are then most consistent with the nominal age of the NGC 6193 cluster, and the anomalously high ${T}_{\mathrm{eff}}=13\,{000}_{-4000}^{+7500}\,{\rm{K}}$ for the secondary can be attributed to its extremely high irradiation from the primary.

Wolk et al. (2008) observed HD 149834 (Source 174 in their study) using Chandra for ∼25 hr, starting at UT 2:37 on 2004 October 25. They identified HD 149834 as a flare source using Bayesian blocks (Scargle 1998) with variability detected at 95% significance; its light curve is plotted in Figure 12. Based on the ephemeris we determined in Section 4.2 (see Table 1), the Chandra observations span orbital phases 0.57–0.80, with a formal uncertainty of less than 0.01. However, the possibility of multiple period aliases makes the uncertainty likely larger.

Zoom In Zoom Out Reset image size

Taken at face value, the ephemeris suggests that the observed flare occurs near phase 0.75, when HD 149834 B is accelerating toward us. In any case, along with the UV excess seen in the SED (Section 4.1 and Figure 7), this flare may be a signature of activity from HD 149834. An actual detection of the secondary eclipse in the X-ray light curve would have implicated HD 149834 B as the source.

Following the subtraction of the best-fit model determined in Section 4.2, we investigate the residuals to determine if any remaining signal is present. By comparing the residuals (black) with the background (gray) in Figure 13, we see that the variation induced by the earthshine events is still present, despite the detrending. This event occurs twice, separated by approximately 16 days.

Zoom In Zoom Out Reset image size

We subjected the residuals to an iterative pre-whitening analysis using the pythia python package wherein we calculate a Lomb–Scargle periodogram of the residuals, identify the frequency with the highest amplitude, and fit a sinusoid to the data to determine the optimal frequency, amplitude, and phase of that signal. The optimized signal is then subtracted, and the signal-to-noise ratio (S/N) of the signal is computed by dividing the amplitude of the optimized sinusoid by the average noise level of the periodogram of the residuals after subtraction of the identified signal within a window of 3 d⁻¹. This process is repeated until an extracted peak is found to have S/N < 3.5 At the end of this process, we simultaneously fit all frequencies, amplitudes, and phases to the original data, and recalculate the S/N according to the new residuals. We retain all independent frequencies with S/N > 4 (Breger et al. 1993). Following Breger et al. (1999) and considering the average uncertainty of the extracted frequencies, we retain any combination frequency with S/N > 3.3.

Following Degroote et al. (2009), after extraction, we filter for close frequencies defined as δ f ≤ 2.5f_R, and keep the frequency that has the higher amplitude. Here, the Rayleigh frequency, f_R ≈ 0.036 d⁻¹, is the inverse of the time-base of the TESS light curve. We identify eight frequencies, listed in Table 2, two of which are combination frequencies or orbital harmonics, within the frequency uncertainties. The full Lomb–Scargle periodogram of the binary-subtracted residuals is shown in Figure 14, where 4 and 3.3 times the average noise level are depicted by the dashed–dotted and dotted gray lines, respectively. Furthermore, we identify f₇ as matching both the 34th orbital harmonic. The residuals after binary subtraction, and optimized sinusoid modeling comprised of the frequencies in Table 2, are shown in Figure 13 in black and red, respectively.

Zoom In Zoom Out Reset image size

Additionally, we note that the mass of the primary is within the range of stars that are expected to exhibit internal gravity waves (IGWs) excited at the boundary of the convective core. Both 2D and 3D hydrodynamical simulations demonstrate that IGWs observationally manifest as a low-frequency excess and/or as resonantly excited g-modes in photometry (Rogers et al. 2013; Bowman et al. 2019; Lecoanet et al. 2019; Edelmann et al. 2019; Horst et al. 2020). Given the observed low-frequency excess in the periodogram of the residuals of HD 149834, we fit periodogram according to Blomme et al. (2011) and Bowman et al. (2019) with:

$$\begin{eqnarray}&&{\alpha }_{\nu }=\displaystyle \frac{{\alpha }_{0}}{1+{\left(\nu /{\nu }_{\mathrm{char}}\right)}^{\gamma }}+W,\end{eqnarray} \tag{ 1 }$$

to determine the parameters of the background profile. Here, α₀ and W are the amplitude at ν = 0 and the white noise, respectively. γ controls the slope of the profile, and ν_char is the characteristic frequency at which the profile has half of its initial value. We find that the low-frequency excess is best described with W = 28.6 μmag, α₀ = 78.8 μmag, γ = 2.34, and ν_char = 3.32 d⁻¹. The best-fit profile is plotted in red in Figure 14.

The kinematics, distance, and age of HD 149834 are consistent with it being a member of NGC 6193. The mass of the primary makes it one of the most massive members of this cluster. HD 149834 A is at the mass limit where a star is predicted to produce a supernova at the end of its life, depending on the rotational and internal chemical mixing history of the star. Notably, HD 149834 is not the most massive system in NGC 6193. There are two massive multiple systems, HD 150135 (a binary; Sota et al. 2014) and HD 150136 (a triple; Mahy et al. 2012; Sana et al. 2012; Mahy et al. 2018), both of which are confirmed luminous X-ray sources (Skinner et al. 2005). Using a combination of radial velocities (Arnal et al. 1988; Mahy et al. 2012; Sana et al. 2013; Cox et al. 2017), photometry, and long time-base interferometry (Sanchez-Bermudez et al. 2013; Le Bouquin et al. 2017), Mahy et al. (2018) found HD 150136 to comprise a ∼42.8+29.5 M_⊙ close inner binary and ∼15.5 M_⊙ wide tertiary. Furthermore, they determine a common age of the components to be 0–2 Myr, which is about half the nominal age estimate of the cluster.

As mentioned in Section 4.3, there is some tension between the metallicity and age of HD 149834, which may be explained either by the anomalously hot secondary or by the system being a member of the cluster halo, which is estimated to be slightly older than the cluster core. This, however, puts more contention between the age estimates of HD 150136 and HD 149834. Alternatively, such an age gradient between halo and core members of the same cluster has been observed in other young clusters and is evidence for outside-in star formation scenarios (Getman et al. 2014).

Of the eight extracted signals, six are independent and two are low-order linear combination frequencies or harmonics. The extracted frequencies occur in two distinct frequency regimes. The single independent high frequency we extracted is consistent with the expected frequency range for p-mode pulsations in β Cep stars. Those lower-frequency peaks we extracted occur in a frequency regime that is consistent with rotational modulation, low-frequency instrumental noise, low-order heat-driven p- and g-modes, as well as IGWs and pulsation modes resonantly excited by IGWs. The 9.7 M_⊙ primary is situated in the mass range where both g- and p-mode pulsations are theoretically expected to be excited (Walczak et al. 2015; Moravveji 2016; Szewczuk & Daszyńska-Daszkiewicz 2017).

Furthermore, the values we determined for γ and ν_char are consistent with the values expected for IGWs in massive stars as determined from simulations (Edelmann et al. 2019; Horst et al. 2020) as well as a large sample of OB-type stars observed by Kepler, CoRoT, and TESS (Bowman et al. 2020). While subsurface convection has been proposed to generate similar signals (Cantiello et al. 2009; Lecoanet et al. 2019), it is expected to produce low-frequency profiles with significantly different values of γ, i.e., γ ≥ 3.5 (Couston et al. 2018), than those expected to be produced by IGWs, i.e., 1 ≤ γ ≤ 3.5 (Edelmann et al. 2019). To this end, we find that the low-frequency power excess observed in the TESS data of HD 149834 to be consistent with the presence of IGWs in the massive primary.

There is much debate concerning the theoretical instability strips of SPB and β Cep pulsators. Several studies have demonstrated that the edges of the instability strips can change drastically depending on the choice of metallicity (Daszyńska-Daszkiewicz et al. 2013), opacity table (Salmon et al. 2012; Walczak et al. 2015; Moravveji 2016), and rotation rate (Szewczuk & Daszyńska-Daszkiewicz 2017). One of the most important implications of this is that these stars require a sufficiently high (nearly solar) metallicity in order to excite pulsations via the κ-mechanism. Furthermore, Southworth et al. (2020) have demonstrated that in addition to the canonical heat-driven pulsations present in β Cep stars, B-type stars in binaries can also host tidally perturbed oscillations. However, we note that due to the low amplitudes expected in such pulsations, we are not likely to observe them with the current data set. To further complicate this picture, both 2D and 3D simulations have demonstrated that IGWs can resonantly excite p- and g-modes (Edelmann et al. 2019; Ratnasingam et al. 2020). With the current frequency resolution, however, we are not able to distinguish between coherently driven modes excited via the κ-mechanism and modes resonantly excited via IGWs.

The evolutionary status of the primary is consistent with the instability regions of both low-order β Cep-type p- and - modes, as well as high-order SPB-like g-mode pulsations. (Godart et al. 2017). Depending on the choice of metallicity, opacity enhancement, and/or rotation rate, p-mode pulsations are only predicted toward the middle to latter half of the main sequence. Based on the evolutionary stage of the primary, the detection of the p-modes is consistent with low-overtone non-radial p-modes. Additionally, both p- and g-modes are expected to be shifted to higher frequencies in the presence of rapid rotation, making the observed g- and p-mode pulsations consistent with low-overtone non-radial oscillations.

Burssens et al. (2020) recently reported on the detection of SPB, β Cep, and hybrid pulsators observed in the TESS 2 minute cadence data. Of the 98 stars in their sample, only nine are eclipsing binaries, and only three are hybrid pulsators. Interestingly, none of the EBs are observed to host hybrid pulsating stars. To date, there are only a handful of hybrid SPB/β Cep stars known in binary systems (see, e.g., Degroote et al. 2012; González et al. 2019) and a handful of either SPB or β Cep pulsators in systems that are eclipsing (see, e.g., Freyhammer et al. 2005; Tkachenko et al. 2014; Jerzykiewicz et al. 2015). While there are several known β Cep or SPB pulsating stars in clusters (Balona & Shobbrook 1983; Handler et al. 2007; Saesen et al. 2013; Moździerski et al. 2019), no hybrid pulsating SPB/β Cep stars in EBs have been identified in clusters to date. Although, this is likely a consequence of data quality and duty cycle, rather than a physical deficiency.

In addition to the pulsational variability, any rotational modulation is also expected to occur in the low-frequency regime. To test this, we consider two scenarios: (i) rotation synchronous with the orbital period, i.e., f_rot = 0.21759 ± 0.00001 d⁻¹, and (ii) rotation consistent with the projected rotational velocity measured by Huang & Gies (2006b), $v\sin i=216\pm 11\,\mathrm{km}\,{{\rm{s}}}^{-1}$, i.e., f_rot = 0.88 ± 0.05 d⁻¹, using the values for R₁ and i from Table 1. Both of these values are plotted in Figure 14, with the rotation rate derived assuming synchronous rotation in red, and the rotation frequency derived from Huang & Gies (2006b) in gray. There is no significant signal remaining in the residuals at the orbital frequency. However, f₄ in Table 2 agrees with the derived rotation rate from Huang & Gies (2006b), within 2σ. Assuming f_rot = 0.88 ± 0.05 d⁻¹, we find v/v_crit = 44 ± 4%, using the values from Table 1. The rotation velocity of the primary is consistent with the observed distribution of projected rotational velocities of B-stars in clusters (Garmany et al. 2015).

The potential identification of a photometric signal associated with the rotation rate raises further questions. Given the stably stratified radiative envelope predicted in a typical 9.7 M_⊙ star, starspots are not expected. Instead, if a surface inhomogeneity is observed, it could be caused by a strong magnetic field that induces surface chemical peculiarity (Alecian et al. 2015; David-Uraz et al. 2019).

An alternative mechanism for inducing photometric variability that coincides with the rotation period is rotationally modulated clumpy winds, as has been observed by Aerts et al. (2018). To frame the context of HD 149834, only about ∼10% of OB stars are observed to have a detectable surface magnetic field. Furthermore, whereas nonmagnetic chemically peculiar late-B-type stars such as HgMn stars are observed to have a high binary incidence (Takeda et al. 2019), higher-mass magnetic stars in close binaries are far less common (Alecian et al. 2015). However, the classification of the primary as magnetic or chemically peculiar requires high signal-to-noise spectra and spectro-polarimetric observations, which are beyond the scope of this paper.