© ESO, 2017

1. Introduction

Massive stars are important for the chemical enrichment of the Universe. Slowly pulsating B (SPB) stars are not amongst the most massive stars, but they share a similar internal structure and are therefore ideal to improve our understanding of massive stars.

Slowly pulsating B stars (SPB) were introduced to the zoo of variable stars by Waelkens (1991). They are non-radial multi-periodic oscillators on the main sequence between spectral type B3 and B9, with an effective temperature ranging from about 11 000 to 22 000 K, and a mass between 2.5 and 8 M_⊙ (e.g. Aerts et al. 2010). They oscillate in high-order gravity (g) modes with frequencies typically ranging from 0.5 to 2 d^-1, which are driven by the κ-mechanism acting due to the iron-group element opacity bump (e.g. Dziembowski et al. 1993). Consecutive radial order n gravity modes of the same spherical degree l are expected to be equally spaced in period, and deviations from this regular pattern carry information about physical processes in the near-core region (e.g. Miglio et al. 2008).

SPB stars are expected to be dominated by a convective core and a radiative envelope, and therefore experience internal mixing processes, which have a significant influence on the lifetime of the star by enhancing the size of the convective region in which mixing of chemical elements occurs. Such a mixing might be induced by convective core overshooting but also by internal differential rotation (e.g. Aerts et al. 2003; Dupret et al. 2004). Despite their importance for realistic stellar structure and evolution models of massive stars, the physical details describing these processes are hardly known. This is mainly because of the small number of detailed investigations of SPB stars, partly due to the few identified modes in these studies (for a recent review see Aerts 2015).

The Canadian space telescope MOST (Walker et al. 2003; Matthews et al. 2004) was very successful in providing the high quality data of SPB stars that are necessary to challenge theory (e.g. Walker et al. 2005; Aerts et al. 2006; Gruber et al. 2012; Jerzykiewicz et al. 2013). The breakthrough in observing SPB stars came, however, with the Kepler mission. Only recently, Pápics et al. (2014, 2015) reported on the detection of a rotationally affected series of g-modes in the two SPB stars KIC 7760680 and KIC 10526294 that show clear signatures of chemical mixing and rotation and which enabled the first actual seismic modelling of SPB stars. A limitation in this respect is that Kepler can only observe fairly faint stars, for which additional observational constraints (e.g., from high resolution spectroscopy or interferometry) are difficult to obtain.

HD 201433 (HR 8094, V389 Cyg) is one of the brightest stars (V ≃ 5.61 mag) suspected to be a SPB star (due to its position in the Hertzsprung-Russel diagram). The B9V star is known to be member of a single-line spectroscopic triple system (Barlow 1989). The most recent determinations of T_eff = 12 193 ± 360 K, log g = 4.24 ± 0.2, and vsini = 15 km s^-1 were published by Takeda et al. (2014). Its Hipparcos parallax is 8.64 ± 0.55 mas (van Leeuwen 2007), from which we obtain an absolute visual magnitude of M_V = 0.29 ± 0.14  mag. Interpolation in the tables of Lejeune & Schaerer (2001) for [Fe/H] = 0.0 (see Sect. 9.3) indicates a bolometric correction of BC_V = −0.70 ± 0.05. With M_bol,⊙ = 4.76 mag (Kopp & Lean 2011) we then obtain L/L_⊙ = 115 ± 15. These parameters locate HD 201433 close to the cool border of the SPB domain.

In this paper we report high-precision photometric observations of HD 201433 with BRITE - Constellation¹, which is an array of five nanosatellites devoted to high-precision, long-term photometry of bright stars as is described by Weiss et al. (2014). Our photometric analysis is primarily based on 156 days of BRITE-Toronto (BTr) observations supplemented by about 13 days of BRITE-Lem (BLb) data (see Table 1). The data products and necessary post-processing are described in Sects. 2 and 3. The Bayesian frequency analysis (Sect. 4) reveals a sequence of nine significant close pairs of frequencies, consistent with rotationally split dipole modes, from which we extracted an average period spacing and rotational splittings (Sect. 5). In Sect. 6 we demonstrate that our interpretation of the BRITE photometry is fully consistent with the signal found in the almost eight-year long SMEI observations. We then construct a dense stellar model grid and search for a representative model of HD 201433 (Sect. 7), which we use in Sect. 8 to infer the internal rotation profile. To complement the space photometry we obtain new high-resolution spectra which extend the time base to slightly more than 96 yr with a total of 231 spectra usable for an orbital analysis. Based on an entirely Bayesian analysis of the radial velocity measurements we improve the published orbital elements and find evidence for a continuously decreasing orbital period. Putting this in context of our asteroseismic and spectroscopic results we conclude that the main component of HD 201433 is a SPB star of about three solar masses showing no significant rotational gradient throughout most of its interior. However, we find indications for tidal interaction with a close companion, causing an acceleration of the outermost envelope. We discuss our findings for HD 201433 in a broader astrophysical context in Sect. 10 and summarise our analysis in Sect. 11.

2. BRITE photometry of HD 201433

The photometric observations used in this study were carried out with two of the five BRITE - Constellation satellites. Each of the 20 × 20 × 20 cm satellites hosts an optical telescope of 3 cm aperture, feeding an uncooled CCD, and is equipped with a single filter. Three nanosats have a red filter (550–700 nm) and two have a blue filter (390–460 nm). The orbital periods are close to 100 min, enabling continuous observations of the chosen target fields for about 5–30 min per orbit.

The detector is a Kodak KAI-11002M CCD with about 11 million 9 × 9 μm pixels (plate scale of 27.3′′ per pixel), a 14-bit A/D converter, an inverse gain of about 3.5 e⁻/ADU, and a readout noise and dark current of about 16 e⁻ and 20 e⁻/s per pixel, respectively, at +20° C. The saturation limit of the pixels at this temperature is about 13 000 ADU, with the response being linear up to about 9000 ADU. Further details about the detector and data acquisition of BRITE - Constellation are described by Pablo et al. (2016)

A problem affecting the BRITE nanosatellites is the higher-than-expected sensitivity of the CCDs to particle radiation, which posed a major threat to the lifetime and effectiveness of the BRITE mission. The impact of high-energy protons causes the emergence of hot and warm pixels at a rate much higher than originally expected. The affected pixels more easily generate thermal electrons and thereby significantly impair the photometric precision of the observations. An additional important problem that appeared after several months of operation was the charge transfer inefficiency also caused by the protons. There were serious problems in the early phase of the mission, but thanks to slowing the readout time and adopting a chopping technique for data acquisition, the effect of CCD radiation damage on the photometry is now significantly reduced (Popowicz et al. 2017).

Satellite pointing is adjusted slightly between consecutive exposures in the chopping mode, so that the target PSFs alternate between two positions (about 20 pixels apart) on the CCD. This means that the PSF-free part of a given subraster image acts as a dark image for the subsequent exposure and subtracting consecutive exposures results in an image with one negative and one positive target PSF. The background defects are thereby almost entirely removed. More details about this technique are given by Popowicz et al. (2017).

HD 201433 was one of the targets in the BRITE - Constellation Cygnus II field and was observed with BRITE-Toronto² for about 156 days in June–November 2015 typically 48 times per BRITE orbit with an average cadence of 5 s exposures every 20.3 s. A significantly shorter data set was obtained with BRITE-Lem³, which observed HD 201433 for about 13 days in September 2015 for typically 30 times per orbit (see Table 1).

As the stellar flux is extracted from differential images (chopping mode) no bias, dark, and background corrections are necessary. The remaining main step is to identify the optimal apertures and to extract the flux within these apertures (Popowicz et al. 2017). The light curves resulting from this pipeline reduction are deposited in the BRITE - Constellation data archive from where we extracted the data of HD 201433 and applied some post-processing routines, as are described in the following.

3. BRITE data post-processing

The raw stellar flux still includes instrumental effects and obvious outliers and therefore needs some post-processing. The BTr data came in two different setups⁴ and five data blocks of approximately equal length, which we treated independently. The first setup at the beginning of the observing run addressed 24 stars in the field but had to be reduced to 18 stars (2nd setup), because of data transfer limitations. Subdivision of the dataset into blocks was required due to a limit of typically 30 000 frames for the standard data reduction software.

Adapting the recipe of Pigulski et al. (2016) we perform the following steps for each of the five blocks:

• Divide the data set into two sets corresponding to the alternatingPSF position on the CCD subraster. (see Fig. 1).

• Compute a 2D histogram of the X/Y positions and fit a 3D multivariate Gaussian to it. Measurements that were obtained with the PSF centre positioned outside three times the widths of the Gaussian (see Fig. 1) are eliminated from further processing. This procedure identifies most of the outliers (about 2.5% of the original data) and rejects them.

• The procedure resumes with the “cleaned” data set and applies a 4σ-clipping to the whole data set, where σ was determined from the complete set. The procedure results in the elimination of additional ~0.3% of all data points.

• To better access the instrumental correlations we first pre-whiten the two highest amplitude frequencies (see Sect. 4), which are subsequently added back after post-processing of the data.

• In the case of HD 201433 the instrumental flux increases typically by 5–7 ADU/s per °C with increasing CCD temperature. A quadratic fit with the CCD temperature is sufficient to correct for this temperature correlation (see top panel of Fig. 2).

• Pixel-to-pixel sensitivity variations of the detector are reflected in correlations between the instrumental flux and the PSF position on the CCD (about 3–5 ADU/s per pixel). We correct for this with polynomial fits (middle panels of Fig. 2).

• Residual instrumental signal is apparent when phasing the instrumental flux with the satellite’s orbital period. We correct for the high-overtone signal with a 200 point box-car filter in the phase plot (bottom panels of Fig. 2).

• The residual instrumental flux is then divided by its average value for conversion to relative flux.

The ten reduced data sets (two for each setup) are then simply stitched together, where no significant offsets at the subset interfaces are found. The final light curve of HD 201433 consists of about 102 000 individual measurements and is shown in Fig. 3. The post-processing reduces the point-to-point scatter of the BTr data of HD 201433 from about 46 to 34 ADU/s (or ~0.9%).

The intrinsic variability of HD 201433 acts on time scales of no shorter than a few hours; hence, the average cadence of about 20.3 s (which corresponds to a Nyquist frequency of ~2130 d^-1) is unnecessarily short. Given this and because the frequency analysis (see Sect. 4) requires good estimates for the uncertainties of the individual measurements we bin the light curve. To keep the dominant cadence short enough (i.e., the Nyquist frequency high enough) we bin the typically 48 measurements per BRITE orbit into two bins, where the standard deviation of the original measurements within a given bin provides a good estimate for the photometric accuracy. The binned light curve consists of about 4200 data points with a median cadence of ~8.13 min (f_nyq ≃ 89 d^-1) and an average error of about 2.1 ppt.

The BLb raw dataset of HD 201433 has a considerably shorter time base than the BTr observations and is due to the higher CCD temperature also much noisier. The original ~1700 measurements reduce to about 1300 useful data points in the post-processed light curve (see Fig. 4). Binning of the typically 30 measurements per BRITE orbit results in 87 data points with a median cadence of ~5.4 min and an average error of about 16 ppt.

4. Frequency analysis

SPB stars are expected to show long-period g modes in a frequency range of up to a few cycles per day. This is well separated from residual instrumental signal and alias peaks due to the orbital frequency of BTr of ~14.7 d^-1 (and multiples of it). We compute the Fourier amplitude spectrum of the unbinned light curve and find no significant peak between 2 d^-1and the Nyquist frequency which cannot be attributed to the satellite’s orbital frequency. The pulsation spectrum (see Fig. 5) of HD 201433 shows the strongest peaks between about 0.5–1 d^-1, and another group of peaks between about 1.5–2 d^-1. Above that no significant power can be found.

An unusual feature in the Fourier spectrum of HD 201433 is that some peaks appear slightly broader than expected from the spectral window function (see inset in Fig. 5). This indicates the presence of close frequencies that are separated by less than (or close to) the formal frequency (Rayleigh-) resolution of 1 /T ≃ 0.0064 d^-1. Such features are problematic for a standard frequency analysis (based on a strict pre-whitening procedure) because assuming a mono-periodic signal in the vicinity of the considered Fourier peak yields a frequency that corresponds to the weighted average of the intrinsic frequency multiplet and pre-whitening this “wrong” signal causes artificial peaks in the spectrum. Furthermore, it is difficult (or often impossible) to objectively rate the significance of the result and its uncertainties.

We use a probabilistic approach to tackle this problem. Kallinger & Weiss (2016) have developed a fully automated Bayesian algorithm that searches for close frequencies in time series data and tests their statistical significance by comparison to a fit with constant (i.e., no periodic) signal and a fit with a mono-periodic signal. The procedure performs the following steps:

• Compute the Fourier amplitude spectrum up to2 d^-1 and determine the frequencywith the highest amplitude.

• Fit N functions, , to the time series, where n incrementally increases from 1 to N so that, in total, N models with 1, 2,..., N sinusoidal components are fit to the data. A, f, and Φ are the amplitude, frequency, and phase of the ith component, respectively. The parameter c serves as an offset to ensure that ^∫_TF_(t)dt = 0 even if the duration T of the time series is not an integer multiple of the signal period. For the fit we use a Bayesian nested sampling algorithm (MultiNest; Feroz et al. 2009), and allow the individual frequencies to vary around the initial frequency by ± 2 /T, and the amplitudes between 0 and 50 times the initial amplitude from the amplitude spectrum. Phases have no initial constraints and can vary from 0 to 1.

• To rate if a signal is statistically significant (i.e., not due to noise) and if so, which model best represents the data, we compute the model probability (p_n) by comparing the global evidences⁵ (z_n) of the fits to those of a fit with a constant factor (z_c). If p = ∑ z_n/ (z_c + ∑ z_n) > 0.95 we consider the solution as real⁶ and not to be due to noise. If so, the best-fit model is then the model with p_n = z_n/ ∑ z_n> 0.95. This means that in order to be accepted, a multiperiodic solution needs to fit the data considerably better than the monochromatic solution. Our approach for the statistical significance of a signal compares well to classical approaches like a SNR> 4 (e.g. Breger et al. 1993; Kuschnig et al. 1997) but has the advantage of providing an actual statistical statement that is based only on the data and that allows us to discriminate between mono- and multi-periodic solutions for closely separated frequencies. In the present case we tested models with up to three components but find that for none of the identified multiplets is a solution with N = 3 statistically significant.

• The best-fit parameters and their 1σ uncertainties are then computed from the marginalised posterior distribution functions as delivered by MultiNest.

• The best-fit model is subtracted from the time series and the procedure starts from the beginning.

We stop the procedure when p drops below 0.66 (corresponding to weak evidence) but we accept only those frequencies with p> 0.95. We note that the frequency, amplitude, and phase uncertainties that are computed from the posterior probability distributions compare well with uncertainties determined from other criteria (e.g. Kallinger et al. 2008).

Based on extensive tests with synthetic data (with the sampling and noise characteristics of the BTr data of HD 201433) Kallinger & Weiss (2016) have shown that the algorithm is capable to reliably (>99.9%) distinguish between a single frequency and a pair of close frequencies if the frequencies are separated by more than ~ 0.5 /T and their amplitudes are larger than about 1 ppt. The uncertainties of the individual frequencies are thereby only slightly larger than for an unperturbed mono-periodic signal but rarely exceed 0.1 /T.

4.1. Frequencies and frequency combinations in the BTr data

Our Bayesian frequency analysis algorithm identified 9 “features” in the binned BTr data of HD 201433 that consist of statistically significant closely separated frequencies in addition to a further 11 single frequencies. An example for a pair of close frequencies is illustrated in the left inset in Fig. 5, where we show the posterior parameter distributions of the one-frequency and two-frequency model fits for the highest-amplitude peak in the BTr spectrum of HD 201433. The evidence of the two-frequency model is orders of magnitude better than for the one-frequency model (despite the Bayesian “penalty” for introducing additional free model parameters), which indicates – based on solid statistical grounds – that more than one frequency is needed to reproduce the data in this frequency range.

The 29 significant frequencies detected in the BTr data set are listed in Table 2. After pre-whitening them from the data, the residual spectrum (see Fig. 5) has an average amplitude of about 110 ppm. The “bump” around 1 d^-1, however, indicates that there is still some undetected signal left (which can well be of instrumental origin). We searched for linear combinations among all significant frequencies. Out of the 29 detected frequencies we find 22 independent frequencies. The remaining peaks correspond to first-order linear combination frequencies (where f_i = f_j ± f_k). In order to be identified as a linear combination a frequency has to fulfil the criterion and its amplitude must be smaller than the amplitudes of its parent frequencies (f_j and f_k). We also searched for higher-order combinations but found none. A schematic view of the independent and combination frequencies is shown in Fig. 6 indicating that all peaks above 1 d^-1 are actually combination frequencies and that 18 of the 19 peaks between 0.4 and 1 d^-1 are part of a pair of close frequencies fully consistent with rotationally split dipole modes.

4.2. Frequencies in the BLb data

The BLb data set has a much shorter time base and is noisier than the BTr dataset. Consequently, the frequency analysis is more challenging (see the amplitude spectrum in Fig. 4). An independent analysis gives only one significant peak with a frequency of 0.852 ± 0.003 d^-1, which represents a weighted average of f₁ and f₄ of the formally unresolved frequencies in BLb. We can, however, fix the frequency to the values determined for the BTr data and fit only the amplitude and phase. We tried various combinations of the three largest amplitude frequencies in the BTr data (f₁, f₄, f₅, f₁ ∧ f₄, f₁ ∧ f₅, f₄ ∧ f₅, and f₁ ∧ f₄ ∧ f₅) and find that a fit with f₁ and f₄ gives the (by far) the best model evidence. The resulting amplitudes and phases are given in Table 3. The two frequencies have very similar amplitude ratios and phase differences in the two BRITE passbands, indicating that they have the same spherical degree (e.g. Daszyńska-Daszkiewicz 2008).

5. Period spacings and rotational splittings

From our list of significant independent frequencies we can identify 9 rotationally split doublets (ν₁–ν₉ in Fig. 6), which indicates that all modes have the same spherical degree of l = 1. Our mode identification is summarised in Table 4, where we assume symmetric triplets with the unobserved (i.e., unresolved) central (m = 0) component being located half way between the | m | = 1 components. There are obviously some modes missing so that measuring the period spacing is not straightforward. If we define the period spacing as dP_i = 1 /ν_i−1−1 /ν_i and assume that between ν₅ and ν₆ two modes are missing (i.e., dP₆ ≃ [1 /ν₅−1 /ν₆] /3) we find the period spacings presented in Fig. 7 as solution S₁. There is a strong gradient of dP with the period, which is not consistent with theory (e.g. Miglio et al. 2008). A more realistic solution (S₂ in Fig. 7) is found when assuming 3 modes missing between ν₅ and ν₆ and one mode missing between ν₂–ν₃, ν₃–ν₄, and ν₄–ν₅. The resulting period spacings are quite uniform with only small deviations from the median value of about 5030 s, which is a typical value for a star like HD 201433 (see KIC 10526294; Pápics et al. 2014). Our assumption of filling missing modes might appear unrealistic but in fact alternating high- and low-amplitude modes were already observed in KIC 10526294, so that it is not surprising that some modes between detected modes fall below the detection limit of our observations (~0.5 ppt). In Fig. 6 we show predictions of the mode sequences based on both solution but note that we cannot distinguish between them based on the observations alone.

The average rotational splittings are determined as ⟨ δf_rot ⟩ = (ν_{+ 1}−ν_-1)/2 (with the indices of ν indicating the m value of the mode). This means that we again assume symmetric triplets, which contradicts the findings of Pápics et al. (2014), but given the missing m = 0 components this is the only possibility we have. Our splittings are, however, similar to what was reported for KIC 10526294 and have an average value of Δf = 0.0028 d^-1. Even the trend of increasing splittings towards higher periods can be observed (bottom panel in Fig. 7), which implies a non-rigid internal rotation profile of the star. For solid body rotation the rotational splitting of l = 1 gravity modes is equal to half the rotation rate of the star, so that we can estimate the average rotation period of HD 201433 to be about 177 d. Since the assumption of a solidly rotating star is likely wrong, this value represents the average rotation period dominated by the near-core region, where the g modes have the largest contribution to the splittings.

6. Coriolis/SMEI photometry of HD 201433

Even though the frequency doublets identified in the BTr data set represent a statistically solid result, as is demonstrated in Sect. 4, the separations of the individual components are close to the formal frequency resolution of the time series. Such a result might be considered at first glance questionable as it disagrees with the Loumos & Deeming (1978) criterion, which requires two close peaks to have a minimal separation of ~1.5 times the frequency resolution to avoid influence on their apparent frequencies when applying the classical Fourier technique. In our case we have the fortunate situation that we can test our claims with a data set long enough to provide the required frequency resolution. These data are provided by the Solar Mass Ejection Imager (SMEI) on board the Coriolis satellite (Eyles et al. 2003; Jackson et al. 2004).

SMEI comprises three wide-field cameras, which are aligned such that the total field of view is a 180deg and about 3deg wide arc, so that a near-complete image of the sky is obtained after about every 102 min orbit. A detailed description of the data analysis pipeline used to extract light curves from these data is provided by Hick et al. (2007). The data of HD 201433 were taken from the SMEI website⁷ and contain about 33 400 measurements covering almost eight years of near-continuous observations. The classical frequency resolution therefore is ten times better than the frequency splittings which we are discussing for HD 201433 (see Table 4). Details about the data set are given in Table 1.

The SMEI photometric time series are subject to strong instrumental effects, such as large yearly flux fluctuations, which obviously are due to an insufficient background correction. As a remedy we phased the data with a period of one year and computed the median flux in 200 phase sub-intervals and find an annual amplitude of about 0.8 mag. We applied an Akima spline fit and computed the residuals to the smooth one-year variation. Back in the time domain one sees outliers, jumps, irregular intensity changes and low frequency variations in the residuals which have to be removed in order to make the very low amplitude frequencies in question detectable. This “cleaning" was achieved via iterative spline interpolation anchored on the median flux in 2–3 day intervals (in other words applying a high–pass filter) by 3–5σ clipping to remove outliers and repeating this procedure 20 times. As a result of this rather arbitrary procedure any signal is gradually suppressed towards low frequencies.

The resulting SMEI time series of HD 201433 is homogeneously sampled without significant aliasing. The median cadence of about 101.6 min results in a Nyquist frequency of ~7 d^-1, which is high enough to cover the intrinsic variability of HD 201433. The Fourier amplitude spectrum of the SMEI data set is shown in Fig. 8. Compared to the BTr spectrum, the formal frequency resolution is more than 18 times better (~0.00035 d^-1compared to 0.0064 d^-1 for the BTr data), while the noise level in the Fourier domain is obviously much higher (~0.53 ppt compared to 0.11 ppt in the BTr data). The photometric precision is, however, sufficient to clearly detect the three largest-amplitude doublets ν₆, ν₇, and ν₉ found in the BTr observations. In fact, the doublets turn out to be triplets with the central component having a smaller (or comparable) amplitude than the wing components. This clearly indicates that the BTr observations are not long enough to resolve all three components of the intrinsic triplets even when using our present frequency analysis method. The frequencies (and even the amplitudes) of the wing components (red vertical lines in the top panels of Fig. 8) agree well, however, with the signal found in the SMEI spectrum. Also the central components, which were estimated from the midpoints of the (BTr) wing components, are in good agreement.

To quantify the agreement we tried to extract the individual frequencies from the SMEI observations but find a classical pre-whitening sequence to be insufficient as the individual frequencies split up in many components indicating strong amplitude modulations. This is already visible in the actual spectrum showing multiple side-lobes around the expected positions of the frequencies. Such a structure reminds of the pattern of intrinsically damped and stochastically excited modes (so-called solar-like oscillations) produce in the Fourier spectrum. This is why we fit a sequence of three Lorentzian triplets to the spectrum, (1)where P is the Fourier power, f_n is the frequency of the central component of the nth triplet and δf_rot,n is its rotational splitting. For simplicity we use a single lifetime parameter τ for all modes. The mode height h_n scales with a geometric factor ζ_m(i), which depends on m and the inclination i at which the pulsation axis is seen. According to Gizon & Solanki (2003) this geometric factor is cos²i for the m = 0 component and 0.5sin²i for the | m | = 1 components of a dipole mode. More details about the detection of Lorentzian profiles are given by, e.g., Gruberbauer et al. (2009). We again use MultiNest for the fit and find a best-fit inclination and mode lifetime of 68 ± 5° and 680 ± 110 d, respectively. The central mode frequencies and rotational splittings are listed in Table 4 and the best-fit sequence of Lorentzian triplets is shown in the middle panels of Fig. 8. The triplets at about 0.88 and 0.98 d^-1 are well defined so that we can test the assumption of symmetric splittings for them. A fit with a modified Eq. (1) (where we allow for individual splittings) indeed shows that the wing components are equally separated from the central component within the uncertainties of about 5%.

We emphasise here that it does not necessarily mean that modes have indeed a stochastic nature, if Lorentzians work well in extracting the mode frequencies from the SMEI observations. In fact, this is very difficult to prove, because one has to show that the signal phase is not coherent, which requires continuous observations that cover many lifetime cycles of the mode. The SMEI data cover about five lifetime cycles (if the modes were stochastic), which is likely not enough to verify a stochastic nature. We do, however, find a strong amplitude modulation for the individual frequencies. To quantify this we fit a sine function to a 720 d long subset of the time series, where we fix the frequency to the value determined earlier by the Lorentzian fit. Moving the subset across the SMEI data in steps of 100 d gives the amplitude (and phase) of a given frequency as a function of time. This is shown in the bottom panel of Fig. 8. We tested various window sizes but always find the same general behaviour. The strongest modulation is found for the triplet at about 0.84 d^-1. While the amplitude of the central component is fairly constant, the wing components vary in amplitude by a factor of more than four. Also interesting is that the variations are periodic with a timescale of about 1500 d and that they are in anti phase. The amplitude minimum of the m = −1 component coincides approximately in time with the largest amplitude of the m = + 1 component, and vice versa. A similar but less pronounced behaviour can be found for the triplet at about 0.98 d^-1. Only the triplet at 0.88 d^-1 is different. The timescale of its modulation is much longer and the wing components vary approximately in phase. The physical origin for this phenomenon is unknown to us but we note that something similar is seen in the Kepler observations of KIC 10526294 (Pápics et al. 2014). Even though the authors did not follow up on this, many modes illustrated in their Fig. 9 show a multiple peak structure typical for amplitude modulation.

Apart from the three triplets shown in Fig. 8 and a single sharp peak at 0.40013 d^-1 (for which we have no explanation at this point) we do not find any further significant variability in the SMEI data. The data set does, however, allow us to verify several assumptions made during the analysis of the BTr observations:

• Our interpretation of the nine pairs of close frequencies in theBTr data as symmetrically split doublets isconfirmed by the independentSMEI observations. Even though we can onlyverify the three largest-amplitude multiplets, the remainingdoublets follow the same statistical criteria and only theiramplitudes are smaller, but still significant in the BTr time series.

• The frequencies of the central components and rotational splittings agree on average within ~1.8σ and 0.8σ, respectively, between BTr and SMEI.

• The interpretation of individual frequencies extracted from the BTr observations as independent oscillations requires approximately stable signal amplitudes. Even though we definitely find amplitude modulations, their timescales are long enough to consider the oscillations in first approximation to be stable during the 156 d long BTr observations. Even for shorter lengths of the subsets, which are used for the bottom panel of Fig. 8, we do not find evidence for significant amplitude modulations shorter than those mentioned above.

• The amplitude modulations also provide a reasonable explanation for the missing modes (Sect. 5) as it might well be that their amplitudes were below the detection threshold of the BTr observations.

Finally we note that we cannot straightforwardly constrain the inclination angle of HD 201433 since the observed frequencies are heat-driven modes. In this case, the 2l + 1 components of a rotationally split non-radial mode are not excited to the same amplitude, contrary to solar-like oscillators (e.g. Gizon & Solanki 2003). However, the formally best-fit value of 68 ± 5°gives an amplitude ratio between the central and wing components which is at least not inconsistent with the observed ones. It seems to be plausible that we see the star more equator-on than pole-on.

7. Asteroseismic analysis

To interpret the observed rotational splittings in terms of internal differential rotation we need a representative stellar model for HD 201433. We therefore compare the observed g modes to non-adiabatic pulsation modes computed with the GYRE stellar oscillation code (Townsend & Teitler 2013) for a grid of non-rotating equilibrium stellar models along stellar evolutionary tracks that pass through the spectroscopic error box (see Table 7 and Sect. 1). The models are calculated with the MESA stellar structure and evolution code (Paxton et al. 2011, 2013). As we are for the time being only interested in a representative model we restrict the models to a single initial chemical composition of (Y, Z) = (0.28, 0.02) and turn off convective core overshooting. We evolve a set of zero-age main-sequence models with masses ranging from 2.9 to 3.225 M_⊙ (with steps of 0.025 M_⊙) until their core hydrogen mass fraction drops below X_c = 0.3. To achieve sufficient resolution along the tracks we limit the evolutionary time steps to 0.2 Myr, which results in about 5200 models with a typical resolution in X_c of 0.002 close to the ZAMS to 0.004 for the most evolved models.

We then compute l = 1 modes with non-adiabatic frequencies ranging from 0.3 to 1.2 d^-1. To search for a best-fit model we compare the observed frequencies (ν_obs) to the theoretical ones (ν_model) by computing the reduced χ² value (e.g. Pamyatnykh et al. 1998; Guenther & Brown 2004), (2)where N and σ_obs are the total number of observed modes and their frequency uncertainties. The typical numerical uncertainty of the model frequencies (σ_model) is estimated from following the frequency of a specific mode (i.e., with a given radial order) during stellar evolution. We thereby assume the actual numerical uncertainty to be of the order of the point-to-point scatter after subtracting a running average. We find a value of about 0.00003 d^-1, which is comparable to σ_obs in some cases and therefore not negligible.

The resulting multi-dimensional χ² space is shown in Fig. 9 and illustrates that while the grid resolution along the evolutionary tracks is sufficient (indicated by the smooth decrease and increase of χ² around the best-fit model when plotted, e.g., as a function of R) the resolution in mass is still too low to find a model whose frequencies agree within the observational uncertainties. We can, however, identify a model that fulfils our requirements (approximately representing the internal structure of HD 201433). The best-fit model has a χ² value of about 6.3, which means that its frequencies matches the observed ones on average within 2.5 (i.e., square root of 6.3) times the average observational errors. The best-fit model parameters are listed in Table 5 and its frequencies and period spacings are compared to the observational values in Figs. 6 and 7, respectively. We further note that increasing the mass resolution would not allow us to improve the situation without simultaneously computing models with various chemical compositions and convective core overshoot parameters. Pápics et al. (2014) and Moravveji et al. (2015) found strong correlations between the mass, the overshoot parameter, and the chemical composition in their seismic analysis of KIC 10526294 (a star that is very similar to HD 201433 in mass and chemical composition, but less evolved), from which we can estimate that turning on core overshooting could potentially increase the mass of the best-fit model by up to 0.1 M_⊙. They do, however, also find that models with a low overshoot parameter fits the observations best. Also the fact that our best-fit model is outside the expected range from spectroscopy in the HR diagram (see Fig. 10) is not really troubling since changing the chemical composition would again slightly change the mass and therefore the position in the HR diagram.

A disadvantage of the χ² method is its inability to provide uncertainties and therefore to set a limit on which models (within the grid) represent the observations and which do not. The “second-best”-fit model M2 (with its position in the HR diagram significantly different from M1) has a χ² of about 7.85 and therefore fits the observed frequencies marginally less well. We therefore follow the approach of Kallinger et al. (2010) to determine the Bayesian model probability and find that for our model grid about 99.9% of the total probability is concentrated in the close vicinity (about ±15 K) of the best-fit model M1 (M1 itself has a p of about 0.15) and that there is only a marginal probability that M2 provides the best representation of the observations (see Table 5). Even though this would already rule out M2 we further investigate it in order to check if a slightly different internal structure (M2 is more evolved than M1 and has a smaller relative core size) affects the further analysis.

8. Internal rotation profile

Building on the measurements of Pápics et al. (2014) for KIC 10526294, Triana et al. (2015) provided the first internal rotation profile of an unevolved intermediate-mass B-type star. They found the star to rotate near its core-envelope boundary with a period of about 71 d and while their seismic data point towards a counter-rotating profile within the radiative envelope they cannot rule out rigid rotation. Such results are key to tackle one of the big open questions in stellar evolutionary theory: the transport of angular momentum inside stars. With the nine rotationally split g modes identified in HD 201433 we are able to provide another example of an internal rotation profile for a star very similar to KIC 10526294.

If we assume that the cyclic rotation frequency Ω depends only on the radial coordinate r, then the frequency splitting δf_rot of a mode with degree l and radial order n can be written as, (3)where I_n,l is the mode inertia and K_n,l(r) gives the unimodular mode kernel, that is a function of the mode’s displacement amplitudes (e.g. Cox 1980). R_∗ is the radius of the model.

Normalised integrated versions of the rotation kernels of HD 201433 are given in Fig. 11 and basically show how a rotationally split frequency is accumulated throughout the star. Obviously, different modes are more or less sensitive to rotation in different regions of the star. The mode ν₉, e.g., gains about 10% of its frequency splitting from rotation in the thin layer above the convective core, where N² spikes (see Fig. 12). The mode ν₇, on the other hand accumulates almost three times more of its rotational splitting in the same region. This “differential” sensitivity to different parts of the star allow us to resolve the stellar rotation profile.

From the fact that the modes ν₇ and ν₉ have practically the same observed rotational splitting we can already conclude that HD 201433 is either a rigid rotator or it rotates much faster in the outer layers than the inner ones (because otherwise the split frequencies of ν₇ and ν₉ would be different). Several inversion techniques have been developed in the past aiming to determine the internal rotation profile of the Sun. The inversion of Eq. (3) is, however, a highly ill-conditioned problem that requires, e.g., numerical regularisation. Beck et al. (2014) found that classical approaches like the RLS method (e.g. Christensen-Dalsgaard 1990) or the SOLA technique (e.g. Schou et al. 1998) are not well-suited for stars with only a few observed rotational splittings (like red giant stars, but also SPB stars) and easily become numerically unstable or it is very difficult to evaluate the accuracy and especially the reliability of the result.

8.1. Forward-modelling of the rotation profile

We therefore follow the forward modelling approach developed by one of us (TK in Beck et al. 2014). The algorithm computes synthetic rotational splittings for a parameterised rotation profile and compares them to the observed splittings. The form of the profile is thereby very flexible. One can, e.g., implement a linear piece-wise model or a functional form (like a Gaussian or a polynomial). The profile parameters are again fitted with the Bayesian nested sampling algorithm MultiNest. These parameters are Ω_n in the case of a zonal model with n zones, or the coefficients of a chosen function. The advantage of this approach is again to provide realistic uncertainties and its capability to compare different models and rate which one best represents the observations (e.g., differential rotation vs. rigid rotation). The algorithm was developed to handle the few available rotational splittings in red giants and has proven to give reliable results that are fully consistent with other methods (e.g. Beck et al. 2014, 2017).

The modes that are accessible to our analysis probe the radiative envelope of HD 201433 from the boundary of the convective core (at a radius and mass fraction of ~0.11 and 0.15, respectively) up to about r = 0.98 R_∗ (m> 0.9999 M_∗), above which all frequencies fall above the Lamb frequency (see Fig. 12), which sets the upper limit for g modes to propagate. We can therefore not expect to directly get information about the core rotation rate and the outer ~2% of the star. We tried various models for the rotation profile to fit the observed splittings, ranging from rigid rotation, over linear piece-wise models (with up to 9 zones), polynomial functions, Gaussians, multi-Gaussians, to Lorentzian and error functions. A selection of the resulting rotation profiles are shown in Fig. 13 and more details are given below:

• (a)

  Rigid rotation: assuming rigid rotation we find a rotationperiod of 201 ± 5 d (0.0050 ± 0.0001 d^-1). However, as can be seen in the inset of panel a in Fig. 13, the synthetic splittings that result from integrating the mode kernels, folded with the rotation profile, do not fit the observed splittings at all. In comparison with other fits the model probability of about 10^-7 is extremely low. We can therefore rule out an entirely rigidly rotating radiative envelope for HD 201433.

• (b)

  and (c) Multi-zonal rotation profiles: we test various zonal models with a fixed position of the zone boundary, i.e., “hard-wired” in the model and find the best formal model probability (p = 0.515) for a fit with a two-zone model with the zone boundary set at a radius fraction of 0.9 (panel b in Fig. 13). The inner and outer zones rotate with a period of 314 ± 43 d^-1(0.0033 ± 0.0005 d^-1) and 14 ± 3 d ( d^-1), respectively. We do not find another model, even with more zones, that comes close in model probability. As an example we show a three-zone model in panel c of Fig. 13, where we add an inner zone (0–0.2 R_∗) but find it to rotate with practically the same rotation rate as the middle zone. This indicates that the data do not support differential rotation in the radiative envelope of HD 201433 below a radius fraction of about 0.9. The very low model probability compared to model (b) results from the additional parameter in the fit⁸.

• (d)

  Two-zone profile with variable zone boundary: in a next step we tried to locate the transition region between the inner slow and outer fast rotating zone. In contrast to the fixed location of the zone boundary in the original two-zone model we now leave it as a free parameter in the fit. The result is shown in panel d locating the zone boundary at a radius fraction of r/R_∗ = 0.93 ± 0.02. While the rotation period of the inner zone remains almost the same, the outer zone now rotates with a period of about 5.5 d (0.182 d^-1) significantly faster than for the original two-zone model. However, due to the rotational splittings containing less and less information about rotation when approaching the surface, the uncertainties start to dramatically increase.

• (e)

  Gaussian rotation profile: so far, the analysis points towards a slowly and rigidly rotating zone that contains almost the entire mass of the radiative envelope, topped by a thin and significantly more rapidly rotating surface layer. The assumption of a sudden increase in rotation speed (by a factor of 20 or so) at the zone boundary seems, however, physically not very plausible. We therefore use a Gaussian rotation profile in which it turns out that we get the best results (i.e., the best model probability) when fixing the centre of the Gaussian to the stellar surface. While the inner rotation period of 292 ± 76 d ( d^-1) agrees well with the two-zone models, the surface rotation rate is 0.30 ± 0.21 d^-1 (3.33 d) – larger than those of the two-zone models. Since the latter do, however, only give average rotation rates within the zones, we consider the models as equivalent. Also the width of the Gaussian (it drops to half its maximum at r/R_∗ = 0.96 ± 0.01) is consistent with the zone boundary of the variable two-zone model. From a statistical point of view the probability contrast between the two-zone models and the Gaussian (about 1:2 for model e:b and 5:2 for model e:d) is not enough to prefer any of them.

• (f)

  Error function rotation profile: even though more physical than a two-zone model, the Gaussian model is not as capable of reproducing the observed rotational splittings, indicating that the transition from the rapidly rotating surface to the slowly rotating interior might not be Gaussian. We therefore fit an error function Ω(r) = Ω₀ + 0.5Ω₁(1 + erf( [r−r₁] /k)) to the observed splittings, where r₁ gives the location where Ω has risen to half its maximum value (Ω₁) and k controls the slope of the increase. The advantage of this function is that it can reproduce very steep to shallow transitions from slow to fast rotation. Interestingly, the fit gives the best-fit parameters Ω₁ = 0.26 ± 0.17 d^-1( d), r₁ = 0.97 ± 0.02, and k = 0.06 ± 0.03, which results in a rotation profile that is very similar to the Gaussian profile. Even though the error function profile is more complex than the Gaussian (it has 4 instead of 3 free parameters) the model probability is 0.112 – high enough to be considered a reasonable representation of the data.

Note that we limit the rotation rates to values between − 1 and 1 d^-1 for all our fits, hence allowing for counter-rotation (as indicated by Triana et al. 2015, for KIC 10526294). However, we did not find a single model that includes a statistically significant counter-rotating zone in HD 201433.

The basic result of the rotation analysis is that the measured rotational splittings are inconsistent with the entire radiative envelope rotating at a constant rate, providing instead strong evidence (in a statistical sense) for a slowly and rigidly rotating envelope topped by a thin and significantly more rapidly rotating surface layer. There is some weak evidence that this layer reaches about 4% down in radius (< 0.001% in mass) and that the transition between the slow and rapidly rotating zone is Gaussian-like. It is interesting to mention that the short rotation period is compatible – within the large uncertainties – with the orbital period of the innermost companion of HD 201433, suggesting angular momentum transfer from the companion star to the outer envelope of HD 201433, causing its acceleration.

An important question is how these results depend on the mode kernels used and therefore the specific best-fit stellar model. We repeated the entire rotation analysis with model M2 and find no significant difference to our original analysis with model M1. The only difference is that the rapidly rotating outer layer seems to reach further down (see panels d to f in Fig. 13), but the deviations are still within the uncertainties.

Finally we test the reliability of our approach, i.e., how well can we reconstruct the internal rotation profile? We therefore compute synthetic rotational splittings based on the M1 mode kernels and arbitrarily assumed rotation profiles (ranging from solid-body rotation, multi-zonal profile, the Gaussians with the centre set to the core and to envelope) and add noise to them. As an example, we show in Fig. 14 the result of our forward-modelling approach for a simulated Gaussian rotation profile, similar to what we find for HD 201433. The simulation demonstrates that we can indeed reconstruct the input rotation profile within the limits of the chosen model. This cannot, however, be directly compared to the real data. In the case of the simulation we know the functional form of the profile and therefore “only” reconstruct its parameters. For the real data, this is not the case and the relatively large uncertainties (compare panels e and f in Fig. 13 to the respective profiles in Fig. 14) do also reflect that the model does not perfectly match the functional form of the real rotation profile. The data do not, however, allow one to better constrain the real form of the rotation profile without over-fitting.

9. Multiplicity of HD 201433

HD 201433 has been known for almost a century to be a gravitationally bound multiple system. Young (1921) measured the radial velocity (RV) of the system 54 times during 1920/21 and found the RV’s to be consistent with a 3.3137 d circular binary orbit. The residuals to the fit were, however, much larger than the uncertainties of the measurements, which is why the system was found to be “peculiar”. Guthnik (1938, 1939, 1942) picked up on the peculiarity and improved the original orbital elements of Young (based on 90 new measurements between 1936 and 1940) but also had no explanation for the persistent large residuals. He also detected photometric variability that was not compatible with the binary orbit. Two Cepheid-like frequencies, 0.838026 d^-1 and 0.885645 d^-1, which are consistent with f₁ and f₄ of the BTr observations (see Table 2), were found in these old photoelectric data, interrupted by periods of irregular variations or constant brightness. Hoffleit (1977) noted that the difference between these two periods is exactly 1/21 d but had no interpretation for this. In fact, 1/21 d corresponds to about 4100 s and is therefore a reasonable estimate for the g-mode period spacing in HD 201433.

Gieseking & Seggewiss (1978) obtained 60 new RV measurements in 1975–77 and improved the original orbital period to 3.3131168 ± 0.000008 d using the combined data set, which was finally large enough to discover the second companion with a period of 154.09 ± 0.02 d, assuming circular orbits. Barlow (1989) re-reduced the observations of Young (1921) and used all available data to find a new solution for the triple system with periods of 3.31311718 ± 0.0000045 and 154.072 ± 0.018 d, and also found that the long-period (LP) companion is in a non-circular orbit with an eccentricity of 0.311 ± 0.057.

9.1. New spectroscopy of HD 201433

Motivated by the fact that the earliest RV measurements date back almost a century and that no high-resolution spectrum of HD 201433 was available for a spectroscopic analysis we organised new spectroscopic observations. A total of 26 spectra (see Table 6) were obtained in May and August 2016 with the Hermes spectrograph (Raskin et al. 2011; Raskin 2011), mounted on the 1.2 m Mercator telescope on La Palma, Canary Islands, Spain. The Hermes spectra cover a wavelength range between 375 and 900 nm with a spectral resolving power of R ≃ 85 000. The wavelength reference was obtained from emission spectra of thorium-argon-neon reference frames in close proximity to the individual exposures. The extraction and the data reduction of the observed stellar spectra were performed with the instrument-specific pipeline (Raskin et al. 2011). Additionally we obtained a single spectrum in late April 2016 with the echelle spectrograph NES (R ≃ 43 000) mounted to the 6 m SAO telescope of the Russian Academy of Sciences.

The radial velocities were extracted by cross-correlating the individual spectra with a template spectrum and fitting a Gaussian to the response function. The RV measurements include the combined orbital motions as well as a component from the surface oscillations of HD 201433 (see Fig. 15). As the latter will influence our orbital solution we have to correct for them. In a first step we increased the formal estimates of the uncertainties by a factor of five ( km s^-1) to make them compatible with the scatter of the first three RV values, which were obtained within a few minutes with Hermes at Mercator and which should practically have the same orbital velocity. The resulting RV uncertainties are still much smaller than the typical uncertainties of the old measurements of about 2.8 km s^-1. The new data should therefore help to significantly improve the orbital elements of the triple system and allow us to search for changes in the orbital periods of the system.

9.2. Line profile variations

A closer look at the individual spectra reveals that the line profiles are in many cases asymmetric and that the profile of a given spectral line changes from spectrum to spectrum. This can be due to the spectral signature of an unresolved companion, but – in our case – is more likely due to non-radial oscillations in HD 201433, which can be best tested with a strong and unblended spectral line. Even though HD 201433 is a relatively hot star, such lines are rare. A good candidate is the Fe ii line at 4233.1 Å. In Fig. 15 we plot the individual spectra corrected for the orbital motions and vertically shifted according to orbital phase (bottom left) and pulsation phase (bottom right) with P_SP from Table 8 and f₁ from Table 2, respectively. If the line profile variations are due to unresolved spectral lines of the companion star the asymmetries in the spectra would be in phase and the “motions” across the line profile would show some structure. This is clearly not the case. If we phase-fold the spectra with the period of the strongest photometric frequency (1.19217 d; Table 2), the asymmetries are now much more in phase. For example, the line profiles of the two spectra at a phase of about 0.55 are almost identical even though they were obtained almost 6 days apart. Furthermore, the line profile variations show some systematics, e.g, a “bump” at a phase of about 0.4 moving towards higher RV half a period later (i.e., at a phase of 0.9). Unfortunately, the spectra are too noisy and the phase coverage is not complete enough for a detailed pulsation mode identification. The available spectra do, however, indicate that the line profile variations are due to non-radial oscillations and that their morphology is not inconsistent with l = 1 modes (e.g. Aerts et al. 2010).

9.3. Atmospheric parameters and v sin i

It is evident that the shape and depth of spectral lines in a single spectrum are strongly affected by oscillations, impeding the determination of atmospheric parameters and the rotational line broadening. We can, however, minimise this effect by averaging all spectra with the additional benefit of increasing the SNR and wavelength resolution. We therefore oversample the individual Hermes spectra, correct them for the orbital motions (original RV in Table 6), and compute a weighted average, where the weights are given by the SNR of the individual spectra. The resulting average Hermes spectrum has a SNR of at about 700 (at 4700 Å).

The atmospheric parameters of HD 201433 were then derived using the SME (Spectroscopy Made Easy) spectral package (Valenti & Piskunov 1996; Piskunov & Valenti 2016), designed to perform an analysis of stellar spectra using spectral fitting techniques in the LTE (local thermodynamic equilibrium) approximation. We choose five large spectral regions 4000–4200 Å, 4200–4400 Å, 4400–4700 Å, 4700–5200 Å, and 5200–5600 Å for the fitting procedure. They include three hydrogen lines, Hδ, Hγ, and Hβ, which are more sensitive to surface gravity than to effective temperature, as well as a large number of lines from light and Fe-peak elements. Few lines of the elements He i, Si ii, Mg ii that are influenced by NLTE (deviations from the local thermodynamic equilibrium) effects were excluded from the fitting. The atomic parameters for all spectral lines were taken from the third version of the VALD database (Ryabchikova et al. 2015a). SME allows one to automatically derive fundamental parameters of a stellar atmosphere: effective temperature, surface gravity, metallicity, radial and projected rotational velocities, micro- and macroturbulent velocities. These parameters are determined by interpolation in a grid of LLmodels (Shulyak et al. 2004; Tkachenko et al. 2012), which are most suitable for the analysis of A to late B-type stars. SME-derived parameters in BAFGK-stars were checked by a comparison with other spectroscopic or independent (e.g., interferometry) determinations (Ryabchikova et al. 2015b, 2016) and were found to agree within 1–2% in effective temperature and ±0.1 dex in surface gravity. We also run SME with Kurucz’s model grid⁹. The final determinations are given in Table 7 together with the results of Takeda et al. (2014), which were derived from calibrations of Strömgren photometric indices.

We find that the spectroscopically derived T_eff and log g values agree well with those derived from photometry. However, the error for the rotational velocity is rather large, obviously due to pulsation. In Fig. 16 we compare the observed line profiles of a rather strong (4555 Å) and weak (4556 Å) Fe ii line with synthetic profiles. While the latter, only broadened by rotation, does not reproduce the observations very well, a combination of rotation and macro-turbulence broadening gives a significantly better fit. There is of course no classical macro-turbulence in the atmosphere of a main-sequence B-star; hence, for HD 201433 the extra broadening is likely due to pulsation. Such an extra broadening is also observed in the atmospheres of, e.g., δ Sct stars (Mittermayer & Weiss 2003) and rapidly oscillating Ap stars (Ryabchikova et al. 2007; Sachkov et al. 2008) but also for SPB stars (De Cat & Aerts 2002). Following the approach of Murphy et al. (2016) we estimate the pulsational broadening to be of the order of 5.5 km s^-1, which reduces the true projected rotational velocity v_esini to about 8 ± 2 km s^-1.

9.4. A new solution for the triple system

With a total of 231 spectra covering slightly more than 96 yr we try to derive a new orbital solution for HD 201433. In a basic Keplerian model for a triple system the primary’s radial velocity RV at an epoch t can be computed according to, (4)where ϕ is the true anomaly, which depends on t, the eccentricity e, the epoch of the periastron T, and P, the orbital period. γ, K, and ω give the systemic velocity, the semi-amplitude of the visible component, and the longitude of periastron, respectively. For a first fit to the observations we again use MultiNest and find a solution which is consistent with the results of Barlow (1989).

We then compute the Fourier power spectrum of the residual velocities (see bottom panels in Fig. 17). There is significant variability in the residuals, but remarkably different for the old and new data sets. While we find a peak at about 1.3130 d^-1with a SNR of 4.8 in the residuals of the old data, there is essentially no power at this frequency in the 2016 data. Contrary to this, we find the strongest peak at about 0.84 d^-1with a SNR of about 3.7 in the 2016 data, which is not detectable in the old data. While this peak can clearly be attributed to pulsation, we have no explanation for the 1.31 d^-1 periodicity. We note that also a 1 d^-1 alias of the SP orbital frequency can be excluded. To improve the orbital solution, we correct for both periodicities according to: (5)where {A, Φ, φ} = {3.22 km s^-1, 0.7616 d, 0.763} and {4.64 km s^-1, 1.1871 d, 0.584} for the old (mBJD < + 24 000) and new (mBJD > + 37 000) data, respectively. The corrected radial velocities from 2016 are listed in Table 6.

In a next step we repeat the fit of Eq. (4) to the corrected radial velocities. As usual the best-fit parameters and their 1σ uncertainties are determined from the posterior probability distributions delivered by MultiNest. The resulting orbital elements are listed in Table 8 and agree with those of the previous fit, but are typically twice as precise as before and are significantly more precise (up to a factor of ten) compared to the solution of Barlow (1989). Most notably, the period of the short-period (SP) orbit is now precise to about 8 ppm (or 0.034 s). This might appear unrealistically small but one has to keep in mind that the data cover more than 10 000 revolutions. While in previous analyses the SP orbit was assumed to be circular we do actually allow also for elliptic orbits, but find the eccentricity to be close to zero (e_SP = 0.015 ± 0.003).

Close binary systems are subject to several physical perturbations, like tidal friction, magnetic braking, mass transfer, etc. (see, e.g., Eggleton 2006) that increase/decrease the angular momentum of the system and therefore affect the orbital period. In fact a number of close binaries are known that show period changes of the order of some ten milliseconds per year for orbital periods of a few hours (e.g. Lohr et al. 2012). The orbital period of the SP companion of HD 201433 is much longer, but given the long time base of the observations we might be able to find a period change. We therefore replace the constant orbital period P_SP in Eq. (4) by the following expression: (6)where P_SP gives the orbital period at T_SP and 1/P·dP/ dt is the relative change rate per year. A fit with MultiNest gives a statistically significant period change rate of (3.2 ± 1.0) × 10^-8. The model evidences of the fits with constant and variable P_SP have an odds ratio of about 1:500, clearly favouring the variable-period fit. Both fits are shown in Fig. 17 and their orbital elements are listed in Table 8. We find only P_LP and e_LP to differ by more than 3σ between the two models, with the variable-period model resulting in a slightly longer and more eccentric orbit for the long-period companion. Remarkably, the rms scatter of the residuals is still much larger than the uncertainties of the RV measurements and the residuals do not show any systematic trends. We therefore attribute the excess in scatter to unresolved oscillations, but tidal effects might also play a role (e.g. Sybilski et al. 2013; Zahn 2013).

The determined period change rate indicates that the orbital period of the short-period companion is now about 0.88 s shorter than in 1921, which is a very small value. But one should keep in mind that the SP period is determined with an accuracy of 0.034 s! Even though we find strong evidence against a stable period we cannot rule out other explanations, like the influence of another star being gravitationally bound to the system. In fact, there is a known visual companion located about 3.4 arcsec from the main component. Its spectral type was indicated between A0 and A7 (Meisel 1968; Abt & Cardona 1984) with an orbital period likely longer than 1000 yr. Clarity might be achieved by searching for variability in the γ-velocity of the triple system. If the detected period change was entirely caused by light-time effects due to orbital motions of an undetected forth component, the expected centre-of-mass velocity change would be on the order of 90 cm s^-1. Obviously, this is well below the accuracy of the present observations (see bottom panels of Fig. 17). On the other hand, a 1000 yr-long orbit would not be detectable in such a Fourier spectrum, even if we had cm s^-1 precision for our 96 yr-long observations. We can, however, search for a trend in the residuals of our analysis. We therefore fit the residuals to a constant period RV time series with first- and second-order polynomials and calculate their MultiNest evidence. These evidence, compared to the case of no trends in RV, i.e. a constant RV, turned out to be statistically insignificant and we therefore can exclude such trends with an amplitude larger than about 30 m s^-1.

9.5. Inclinations of the orbital planes

Because we do not find eclipses or tidally induced flux modulations in the light curve¹⁰, which is larger than about 0.3 ppt at the orbital period of the SP companion, it is not possible to directly infer the inclination (i) of the orbital planes from the available RV measurements. Hence we can only determine the projected semi-major axes and the mass functions (see Table 8): (7)where the indices 1–3 indicate the central star and SP and the LP companion, respectively. From seismology we have, however, a good idea about the mass of HD 201433 and can therefore estimate the companion masses M₂ and M₃ as a function of i. From this we can now determine the orbital separations: (8)where q is the mass ratio M₁/M_2,3. Considering all uncertainties and a conservative estimate of M₁ = 3.0 ± 0.2M_⊙ the resulting companion masses and orbital separations are shown in Fig. 18 as a function of the orbital inclination. In this context we can now test various assumptions.

A possible limitation for the inclination comes from the fact that we do not see any sign of eclipses or rotational modulation in the BRITE photometry, e.g., due to a distorted geometry of a contact or semi-detached binary. For a given orbital separation (a), a limit for i to avoid an eclipse is determined by cosi< (R₁ + R₂)/(a), where we assume , with ζ = 0.57 and 0.8 for stars with a mass above and below 1 M_⊙, respectively (e.g. Torres et al. 2010). The resulting maximum inclination can be translated into a minimum M₂ for which the system does not eclipse. Figure 18 shows this limit as a function of the orbital separation in units of R_⊙ and it is obvious that the system can only eclipse for inclinations larger than about 85°.

Another constraint comes from ellipsoidal variations in non-eclipsing close binary systems whose components are distorted by their mutual gravitation. According to Morris (1985), the peak-to-peak amplitude of the primarys’ light variations can be connected to the system parameters as , where and u₁ and τ are the limb-darkening and gravity-darkening coefficients, respectively. For small amplitudes , with and I₂/I₁ being the semi-amplitude of the ellipsoidal variations and the intensity ratio of the two stars, respectively, in the BTr passband. For stars of similar temperature their luminosity ratio translates directly into an intensity ratio in a given passband, but for stars with significantly different temperatures the different spectral energy distributions need to be taken into account as well. For simplicity we only consider the effects of black-body radiation. For main-sequence stars L ∝ M^3.9 and T ∝ M^{1−ζ/ 2} (with ζ = 0.8 and 0.57 for M larger and smaller than 1 M_⊙, respectively). The mass ratio then translates into an intensity ratio according to I₂/I₁ = (M₂/M₁)^3.9E_(λ,T1)/E_(λ,T2), where E_(λ,T) is the black-body radiation for a temperature T at the central wavelength λ of the passband. Instead of the original bolometric limb- and gravity-darkening coefficients used in Morris (1985) we adopt u₁ = 0.39 and τ = 0.45 from Claret & Bloemen (2011), which are computed for the SDSS r′ passband filter (which is similar to the BTr passband). Ellipsoidal variations usually have a double-peak structure in the light curve so that we can expect to find the strongest peak in the amplitude spectrum at twice the orbital frequency of 0.603644 d^-1. We do indeed find a significant frequency at 0.60535 d^-1 (f₁₀ in Table 2), which differs by less than 2σ from twice the orbital frequency. We, however, identify this as the m = 1 component of the doublet ν₃. Even though this interpretation seems more plausible we can most conservatively assume that the signal amplitude of about 0.7 ppt is entirely due to ellipsoidal variations. Using the relation shown in Fig. 18 we find a maximum inclination of 53° that the system can have to produce ellipsoidal light variations with an amplitude of less than 0.7 ppt. This translates into a minimum companion mass and orbital separation of about 0.4 M_⊙ and 14.1 R_⊙, respectively.

At first sight the spectroscopic observations of HD 201433 are dominated by the spectral signature of the primary star with no obvious evidence for spectral lines originating from one of the companion stars. This indicates that the latter are much fainter than the primary star. For a more comprehensive testimony we search for the spectroscopic signature of the inner companion as follows. Its radial velocity is mirrored to the system velocity with an amplitude scaled by the mass ratio of the two stars. We therefore correct the observed spectra for mirrored radial velocities that result from a given mass ratio and average the spectra that were observed at phases around the RV maxima (i.e., ±0.1 in phase, which covers ten spectra). If the companion were bright enough we would find spectral lines at the expected RV. Compared to looking at individual spectra this approach increases the sensitivity by averaging various spectra and also allows us to distinguish between weak lines of the central star and spectral lines originating from the companion. The absence of spectral lines either means that we did not get the mass ratio right or that the companion’s spectral lines are too diluted to stand out of the noise. We tested companion masses from 0.3 to 3 M_⊙ (in steps of 0.1) but did not find any signature of individual spectral lines from the companion in the wavelength range around the Fe ii line, as is illustrated in Fig. 15. This allows us to define an upper limit for the luminosity of the companion star. The averaged spectra (for each companion mass) scatter in normalised flux by about 0.0025 in the range where we expect a spectral line of the companion, which would then need to be at least 0.0075 deep to to be differentiated from the noise by more than 3σ. Ignoring white dwarfs for the time being, the companion is very likely a low-mass main sequence star very close to the ZAMS (given the relatively young age of HD 201433 for a main sequence star). Given this we assume the intrinsic normalised depth of the genuine spectral line of the companion to be at least 0.5. This results in an intensity dilution factor of at least 67 at the given wavelength. As for the ellipsoidal variations, we also need to account for the stars’ different spectral energy distributions to translate the intensity ratio into a mass ratio. Using mass and effective temperature of our representative model (see Table 5) the companion mass needs to be below 0.97 M_⊙ in order to be undetectable in our analysis. According to Fig. 18 this translates into a minimum inclination of about 21° and a maximum orbital separation of about 14.8 R_⊙. The above considerations do not explicitly account for rotational broadening of the companions’ spectral lines but already include some margin for moderate rotation (in fact the considered spectral lines in low-mass stars are much deeper than 0.5). However, even cutting the minimum dilution factor in half (i.e., allowing for more rotational broadening) does only increase the upper mass limit by about 20%.

The above considerations are only valid for main sequence stars but the companion could also well be a white dwarf (WD). They are orders of magnitudes less luminous than HD 201433. It is therefore practically impossible to see the spectral signature of a white dwarf in a composite spectrum with a B-type star. In a statistical sense a WD-companion would have a mass of about 0.6 M_⊙ (e.g. Tremblay et al. 2016) and very little chance to be more massive than one solar mass.

If we assume the outer layers of HD 201433 to rotate with the orbital period of the SP companion (i.e., synchronised rotation), the equatorial rotational velocity would be of the order of v_rot = 2πR₁/P_SP ≃ 39.6 km s^-1. A comparison with the observed v_esini = 6–10 km s^-1 gives an inclination of 9–15°. Based on the above (conservative) considerations we can, however, exclude such a low inclination. In fact, it would require a 1.5–2.5 M_⊙ companion, which we would definitely detect in the spectra. The obvious solutions for this contradiction are that either the rotation axis of HD 201433 is not aligned with the SP orbital axis or its top-layer rotation is not synchronised. Even though very speculative, we remind that the doublet structure of the oscillation modes indicates that we observe the oscillations more equator-on than pole-on (i.e., the inclination of the pulsation axis is “large”). Assuming the pulsation axis to be aligned with the rotation axis would then favour mis-alignment with the SP orbit. On the other hand, our analysis of the internal rotation profile (Sect. 8) indicates that the surface of HD 201433 rotates with a period in the range of 2–11 d, and does therefore not exclude that the surface rotation is somewhat slower than the SP orbital period (i.e., the equatorial rotational velocity is smaller than for synchronised rotation and therefore the inclination is larger than 9–15°). We think the latter explanation is more plausible. Given that the system is quite young, that circularisation time scales in binary systems are shorter than synchronisation time scales (e.g. Zahn 2013) and that the SP orbit is still not fully circularised (e_SP ≠ 0), it is plausible that the gravitational interaction between HD 201433 and its SP companion is just about to finish circularisation. Hence, spin-alignment is still active and the rotation of the main star is accelerating. A better understanding of this scenario would, however, require detailed modelling of the dynamical history of the multiple system, which is beyond the scope of this paper.

10. Discussion

Rotation is a key process of the evolution of stars, that is still not fully understood. Along with magnetic fields, rotation strongly affects the transport of chemical elements and of angular momentum (AM), and therefore the structure and evolution of stars (e.g. Zahn 1992; Maeder & Meynet 2012). Isolated stars are believed to conserve their total AM throughout evolution from the main sequence to high up the giant branch. During this evolution, their cores (mostly) contract while their envelopes (mostly) expand until they are eventually ejected carrying away AM to the stellar environment. If stars would locally conserve AM throughout their life, their cores would spin up and the surviving compact objects would spin much faster than is actually observed (e.g. Kawaler 2015). This implies that standard formulations of stellar structure and evolution miss one (or more) essential processes. Which and what their individual contributions are, as well as over which timescales they operate, is still unclear and none of the proposed mechanisms (e.g. Tayar & Pinsonneault 2013; Cantiello et al. 2014; Fuller et al. 2014; Rüdiger et al. 2015) can sufficiently reproduce the observations.

A first step to improve our understanding of stellar rotation is to provide a clear view on how AM is distributed inside stars as a function of various parameters (e.g., stellar mass) and during different evolutionary stages. Prime candidates for such studies are sub-giant and red-giant stars for which the contrast between the core and surface rotation has been measured using rotationally split mixed p/g modes (e.g. Beck et al. 2012; Mosser et al. 2012). These measurements are, however, only one piece of the puzzle as the star’s rotational history remains unknown and it is thus difficult to compare the observations to model predictions. On the other hand, the outer envelope contains almost all of the star’s AM. Based on the measured surface rotation rate of an evolved star one can thus well approximate how fast its surface was rotating when the star was still on the main sequence by simply following the internal mass distribution of a representative model backwards. By assuming rigid rotation on the MS the entire rotational history of an evolved star can then be normalised to its MS properties, which allows a direct comparison of different stars. We have done this for various stars in Fig. 19, where we show the core-to-surface rotation contrast as a function of the stellar radius for six sub-giants (Deheuvels et al. 2014) and the more evolved red giants KIC 4448777 (Di Mauro et al. 2016) and KIC 9163796 (Beck et al. 2017). Once the diagram is more populated (and such analyses are underway) it will allow us to follow the AM redistribution in the interior of these stars and therefore provide the observational constraints necessary to test the various AM transport mechanisms.

Our concept of “visualising” the AM transport in evolved stars depends critically on the assumption that stars rotate rigidly in first proximity when leaving the MS. This is in fact difficult to test observationally and has so far only been shown for three late A to early F-type stars (KIC 11145123 – Kurtz et al. 2014; KIC 9244992 – Saio et al. 2015; KIC 7661054 – Murphy et al. 2016). Even though Triana et al. (2015) could not rule out rigid rotation for the B8V star KIC 10526294, our analysis of HD 201433 provides the first confirmation that this important assumption also holds true for the progenitors of more massive red giants.

The models used to place the observations in Fig. 19 were computed with the Yale stellar evolution code (YREC; Guenther et al. 1992; Demarque et al. 2008) for near-solar composition and calibrated mixing length parameter (Z = 0.02, Y = 0.27, α_MLT = 1.8) assuming standard solar mixture (Grevesse et al. 1996). More details about the constitutive physics are described by, e.g., Kallinger et al. (2012) and references therein.

11. Summary

HD 201433 has the potential to be a Rosetta-stone SPB star. It is very bright and, hence, allows for detailed investigations with spectroscopic, polarimetric, interferometric and photometric observations on all relevant scales of temporal and photometric resolution, in the visible region and beyond. It is a success of BRITE–Constellation to have identified the potential of this object.

HD 201433 is a known single-lined spectroscopic triple system consisting of a massive SPB star orbited by two low-mass stars with periods of about 3.31 and 154 d. Our 27 new high-resolution spectra, obtained in 2016, extended the time base to slightly more than 96 yr with a total of 231 spectra usable for an orbital analysis. Such a long time base allowed us to check with Bayesian techniques for variability of the short-period companion, and, indeed, we find an annual relative period change rate of (3.2 ± 1.0) × 10^-8, or that the orbital period is now about 0.9 s shorter than about 10 000 orbits earlier, in 1921. Even though we cannot exclude other explanations (like the influence of an undetected third companion) for this it appears plausible that the decreasing period is the signature of an angular momentum transfer in the inner binary system. The orbital elements are presented in Table 8 and agree with Barlow (1989) within the errors, but are now more precise by a factor of 2 to 10. In addition, we also allow for elliptical orbits, but find the eccentricity of the inner companion to be small (e_SP = 0.015 ± 0.003) and close to the circular orbit, as assumed in pervious analyses. Because we do not find eclipses or tidally induced flux modulations in the BRITE observations and since the spectral signature of the companions is not visible in our spectroscopic time series we can limit the mass of the inner companion to the range of about 0.4–1 M_⊙.

Obviously, pulsation effects in spectra interfere with the dynamical analysis of HD 201433 and, hence, a photometric frequency analysis was performed contemporaneously. Indicated by peaks in the Fourier spectrum of the BRITE observations being broader than expected from the spectral window function, we applied a Bayesian-based technique to identify frequencies which are separated by less than the formal (Rayleigh-) frequency resolution of 1/T, with T being in our case 156 d. We can reliably distinguish between a single frequency and a pair of frequencies (Kallinger & Weiss 2016), if the latter are separated by more than ~0.5 /T and their amplitudes are sufficiently large, in our case ~1 ppt. Remarkably, the frequency errors are only slightly larger than for a mono-periodic signal, but rarely exceed 0.1 /T (i.e., one tenth of the Rayleigh frequency resolution).

The pulsation spectrum (see Fig. 5) of HD 201433 shows the strongest peaks between about 0.4 d^-1 and 1 d^-1, and another group of peaks between about 1.5 d^-1 and 2 d^-1. From our list of significant independent frequencies we can identify nine rotationally split doublets (ν₁–ν₉ in Fig. 6), which indicate that all modes have the same spherical degree of l = 1 (Table 4). The available spectra strongly indicate that the line profile variations are due to non-radial oscillations and that their morphology is consistent with l = 1 modes. Assuming that some frequency pairs could not be identified in the BTr data due to insufficient SNR, the average period difference for the sequence of pairs is about 5030 s, which is a typical value for a star like HD 201433. Our interpretation of the BTr photometry is fully consistent with the almost 8-yr long but poorer-quality observations obtained by the SMEI instrument on board the Coriolis satellite, for which we identify three rotationally split triplets whose central frequencies and splittings agree well with those extracted from the BTr data. We also find the amplitudes of the individual components to be modulated on timescales of about 1 500 d.

A highlight from the BRITE photometry of HD 201433 is a trend of increasing frequency splittings towards higher periods, which implies a non-rigid internal rotation profile – an outstanding detection for the group of SPB stars and a key to tackle one of the big open questions in stellar evolutionary theory, the transport of angular momentum inside stars. For this investigation we computed about 5200 MESA models (see Fig.10) and their non-adiabatic pulsation modes with the GYRE code. Among those we search for a representative model that reproduces the observed frequencies best using classical χ² techniques but also more reliable statistical methods. Our best fit model has a mass and radius of 3.05 M_⊙ and 2.6 R_⊙  respectively, and is about half way through its main-sequence lifetime.

HD 201433 is almost a twin star to KIC 10526294 (Pápics et al. 2014). Apart from their age (HD 201433 appears to be more evolved) they are very similar in mass, surface gravity, and effective temperature and even the measured rotational splittings are similar in value and trend. Triana et al. (2015) investigated the internal rotation profile of KIC 10526294 and found the rotation rate near the core-envelope boundary to be well constrained. Their seismic data are consistent with rigid rotation but a profile with counter-rotation within the envelope has a slight statistical advantage over constant rotation. Our analysis of the internal rotation profile is more conclusive and we can statistically rule out counter-rotation in HD 201433.

The modes that are accessible to our analysis probe the radiative envelope of HD 201433 from the boundary of the convective core (r ≃ 0.11 R_∗) up to about 0.98 R_∗ (Fig.12). According to our analysis, the measured rotational splittings do not support differential rotation in the radiative envelope of HD 201433 below a radius fraction of about 0.9, but are, on the other hand, inconsistent with the entire radiative envelope rotating at a constant rate. Our Bayesian analysis provides strong evidence for a slowly (292 ± 76 d) and rigidly rotating envelope, topped by a thin and significantly more rapidly rotating surface layer reaching down to about 96% of the radius. It is interesting to mention that this rotation at the surface is compatible, within the given uncertainties, with the orbital period of the innermost companion of HD 201433, suggesting an angular momentum transfer from the companion to the outer envelope of HD 201433, causing its acceleration.

Combining the asteroseismic inferences with the spectroscopic measurement of the projected rotational velocity and our orbital analysis of the inner binary system we conclude that tidal interactions between the primary SPB star and its inner companion have almost circularised the orbit but not yet aligned all spins of the system and just started to accelerate the rotation of the SPB star to synchronise rotation. This makes the system quite interesting in order to study tidal interaction of binaries “in action”.

We further refer to a detailed follow-up study about the atmosphere of HD 201433. The star’s atypically sharp spectral lines make it a perfect candidate to study NLTE effects in the atmosphere of a hot star (Ryabchikova et al, in preparation).

---

Acknowledgments

We thank the technical team as well as the observers of the Hermes spectrograph and Mercator Telescope, operated on the island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. BTr operations are supported through a Canadian Space Agency (CSA) Academic Development grant. The Polish contribution to the BRITE mission is supported by the NCN grant 2011/01/M/ST9/05914. The research leading to these results has (partially) received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 670519: MAMSIE). T.K. and W.W. are grateful for funding via the Austrian Space Application Programme (ASAP) of the Austrian Research Promotion Agency (FFG) and BMVIT; and W.W. for support of the University of Vienna (IS 538001, IP 538007). P.G.B. acknowledges the ANR (Agence Nationale de la Recherche, France) program IDEE (No. ANR-12-BS05-0008) “Interaction Des Etoiles et des Exoplanetes”. P.G.B. also received funding from the CNES grants at CEA. Y.P. and T.R. thank RFBR grant 15-02-06046 for partial financial support. T.L. acknowledges funding of the Austrian FFG within ASAP11 and by the grant WTZ CZ 15/2017. A. Pigulski and A. Popowicz acknowledge support from the Polish National Science Centre (grant No. 2016/21/B/ST9/01126 and 2016/21/D/ST9/00656). E.S. is thankful to the Russian Science Foundation (grant No. 14-50-00043) for financial support. G.A.W. acknowledges Discovery Grant support from the Natural Sciences and Engineering Research Council (NSERC) of Canada. A.F.J.M. is grateful for financial aid from NSERC (Canada) and FQRNT (Quebec). G.H. acknowledges financial support by the Polish NCN grant 2015/18/A/ST9/00578. K.Z. acknowledges support by the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF, project V431-NBL).

References

All Tables

All Figures