© ESO 2019

1. Introduction

In 1978, the Einstein observatory detected for the first time X-rays associated to massive stars (Harnden et al. 1979; Ku & Chanan 1979; Seward et al. 1979). Among these was Cyg OB2 #12 or Schulte 12 (V = 12.5), a highly reddened and peculiar star known for several decades (Morgan et al. 1954). Schulte 12 was classified as an early B hypergiant of spectral type B3–4Ia+ (Clark et al. 2012). With its extreme brightness, the star would be placed in the Hertzprung–Russell diagram above the Humphreys–Davidson limit. This explains why Schulte 12 was often considered as a luminous blue variable (LBV) candidate (Clark et al. 2005). In support of this classification, evidence for changes in spectral type (B3–B8) have been presented by Massey & Thompson (1991) and Kiminki et al. (2007). Variations of the temperature were later refuted by Clark et al. (2012) because of the low signal-to-noise ratio of the data, prohibiting changes to reach a low significance level, but also because changes in apparent spectral type may not necessarily reflect a varying photosphere hence could agree with a constant temperature. A lack of helium enrichment is also at odds with an evolved star status (Clark et al. 2012). Furthermore, no surrounding ejected nebulosities seem to exist (Kobulnicky et al. 2012; Oskinova et al. 2017).

The early X-ray detection poses a puzzle which could not be solved in the four decades elapsed since then. Indeed, the X-ray emission of Schulte 12 is bright (log[L_X/L_BOL] = − 6.1) and presents a hard (kT  ∼  2 keV) component (Rauw 2011). This is clearly at odds with the properties of LBVs, which generally present very faint (undetected) X-ray emissions (Nazé et al. 2012a). Besides, it is also at odds with typical properties of massive stars. The intrinsic X-rays from massive stars with fast winds (thousands of km s⁻¹) are generally thought to come from embedded wind shocks (e.g. Feldmeier et al. 1997), which generate soft (kT  =  0.3 − 0.6 keV) and modestly bright (log[L_X/L_BOL]  ∼   − 7 for O-stars) X-ray emissions (Oskinova 2005; Nazé 2009; Rauw et al. 2015a). In this context, it is interesting to note that Schulte 12 actually has a very slow wind (ν_∞ = 400 km s⁻¹, Clark et al. 2012 or even 150 km s⁻¹, Klochkova & Chentsov 2004), incompatible with high-temperature plasma – though larger velocities have also been proposed; see Leitherer et al. (1982).

There is however another Galactic LBV that displays bright and hard X-rays: η Carinae. In this case, they arise from the face-on collision between the two winds of the components in this massive binary (see e.g. Hamaguchi et al. 2014, and references therein). Similarly, the LBV HD 5980 in the Small Magellanic Cloud also emits bright and hard X-rays linked to a wind–wind collision (Nazé et al. 2018a, and references therein). A telltale signature of such colliding winds is variability (for a review, see Rauw & Nazé 2016). Indeed, the intrinsic strength of the collision changes with the stellar separation in eccentric binaries and the absorption may also vary as the shocked plasma is alternately seen through the two different winds (Stevens et al. 1992). Because they are linked to the orbital configuration, such changes recur with the orbital period. For Schulte 12, variations of its X-ray emission have been reported several times (Rauw 2011; Yoshida et al. 2011; Cazorla et al. 2014) but no periodicity was found. In view of the apparently monotonic decrease in X-ray fluxes from 2004 to 2013, an alternative scenario where X-rays are produced by fading shocks between slow-moving ejecta and fast winds after an eruption could not be rejected (Cazorla et al. 2014).

The possible presence of a companion to Schulte 12 has been considered in the past (Souza & Lutz 1980) and has subsequently been revived in a few recent studies. First, Caballero-Nieves et al. (2014) reported the detection of a close visual companion at 64 mas with ΔV = 2.3. Maryeva et al. (2016) confirmed its presence and constrained its orbital period to 100–200 yr from its slight change of position. These latter authors also found a third object further out at 1.2″ with Δm = 4.8. However, such distant companions are unlikely to produce an X-ray bright wind–wind collision as plasma density would be too small at such large separations. Until now, X-ray bright wind–wind collisions have been found in massive binaries with periods of a few years (up to a decade in 9 Sgr; Rauw et al. 2016), not decades or centuries. Secondly, small line-profile and radial-velocity (RV) changes were reported by Klochkova & Chentsov (2004) and Chentsov et al. (2013) but they could not confirm or rule out the binarity hypothesis (Chentsov et al. 2013). Only a double-lined spectroscopic binary seemed excluded by the data (see also Clark et al. 2012). Finally, Oskinova et al. (2017) studied a high-resolution X-ray spectrum of Schulte 12 and found rather strong f lines in fir triplets from He-like ions. While the Si XIII triplet is compatible with the absence of depopulation of the upper level of the f line, the Mg XI lines suggest some depopulation, though not due to UV radiation but to collisions: there would thus be dense plasma where X-rays arise. This could be compatible with X-ray generation very close to the photosphere or in a wind–wind collision. The latter scenario is further supported by the (broad and unshifted) profiles of X-ray lines and the lack of very strong absorption. Oskinova et al. (2017) further suggested a late-O type for the optically closest companion, if there was a wind–wind interaction between this star and the LBV (which remains to be demonstrated). However, they interpreted the observed X-ray variability as being random and incompatible with a companion on a long (decades) orbit.

While important information has been gained in recent years, the nature of Schulte 12 is still not fully understood. Time is thus ripe to perform an in-depth study, with a thorough X-ray and optical monitoring (described in Sect. 2). The derived results led to several surprises (Sects. 3 and 4 for X-ray and optical/near-IR findings, respectively). Their possible interpretation is discussed in Sect. 5 and then summarized in Sect. 6.

2. Observations and data reduction

2.1. Optical and near-IR domains

2.1.1. Photometry

Photometry of Schulte 12 was obtained in both the Ic and V filters at the private observatory of one of the current authors (F.C.), situated in Vallirana (near Barcelona, Spain) and equipped with a Newton telescope of 20 cm diameter (with f/4.7 and a German equatorial mount). The camera is a CCD SBIG ST-8XME (KAF 1603ME). The exposures were typically of 180 s and 120 s duration for the V and Ic filters, respectively.

The images were corrected for bias, dark current, and flat-field in the usual way using the data-reduction software Maxim DL v5¹. The photometry was extracted with the FotoDif v3.93 software, using as comparison star SAO 49783 (TYC 3157-195-1). Four stars (TYC 3157-1310-1, TYC 3157-603-1, TYC 3161-1269-1, TYC 3157-463-1) were further measured to check the stability of the solution.

2.1.2. Spectroscopy

Schulte 12 was observed 12 times between 2010 and 2017 with the Hermes spectrograph (Raskin et al. 2011) installed on the 1.2 m Mercator telescope on La Palma (Spain). Given the magnitude of the star, between 2 and 7 consecutive exposures of 1800 s (depending on the weather conditions) were summed up to increase S/N. The data were taken in the high-resolution fibre mode, which has a resolving power of R ∼ 85 000. The spectra cover the 4000–9000 Å wavelength domain but because of the high extinction of the star, the range below 5750 Å was not used for the present analysis. The raw exposures were reduced using the dedicated Hermes pipeline and we worked with the extracted cosmic-removed, merged spectra afterwards.

A set of 19 observations were obtained with the Carmenes (Calar Alto high-Resolution search for M dwarfs with Exoearths with near-IR and optical Echelle Spectrographs; Quirrenbach et al. 2018) instrument at the 3.5 m telescope of the Centro Astronómico Hispano Alemán at Calar Alto Observatory (Spain). These observations were granted via the OPTICON common time allocation process (program IDs 2018A/004 and 2018B/001) and through German open time (program ID F17-3.5-002). Carmenes features a dichroic that shares the light between two separate echelle spectrographs: the first one covering the optical domain from about 5200 to 9600 Å at a spectral resolving power near R = 94 600, and the second covering the near-IR from about 9600 Å to 1.71 μm at R = 80 400. All data were obtained in service mode at a rate of about one observation every 10 d during our main observing campaign in May–August 2018, and the exposure times were 4 × 5 min per observing night. The data were reduced using the Carmenes pipeline (Caballero et al. 2016). The Carmenes spectra taken the same night were combined to improve S/Ns.

Table 1 lists the dates of these spectra, as well as their S/Ns around 7525 Å. Because of the strong extinction, typical S/Ns vary from ∼50 at 5800 Å to ∼300 in the near-IR. To obtain the best diagnostics, we focused on wavelength intervals with strong stellar lines and limited contamination by interstellar or telluric lines. However, small contamination could not always be totally avoided. We therefore performed a correction for telluric lines around Hα, He I λ 7065 Å, and He I λ 10 830 Å, as well as Paschen lines at 8500–8900 Å and 1 μm. This was done within IRAF using the telluric template of Hinkle et al. (2000) in the optical range and a telluric spectrum calculated for Calar Alto with MOLECFIT (Smette et al. 2015) in the near-IR. As a last step, the spectra were normalized over limited wavelength windows using splines of low order.

2.2. X-ray domain

2.2.1. XMM-Newton

The Liège group have obtained ten observations of the Cyg OB2 region since the launch of XMM-Newton. The oldest ones were reported in De Becker et al. (2006), Rauw (2011), Nazé et al. (2012b), and Cazorla et al. (2014). Schulte 12 also appears in six additional XMM-Newton exposures centred on PSR J2032+4127. However, in these exposures, it is far off-axis and only in the field-of-view of the MOS2 detectors. Moreover, stray light from Cyg X-3 contaminates the source (further discussed below). All datasets were reduced with SAS v16.0.0 using calibration files available in Fall 2017 and following the recommendations of the XMM-Newton team². No barycentric correction was applied to event arrival times.

The EPIC observations were taken in the full-frame mode and with the medium filter (to reject optical/UV light), except for ObsID 0677980601 where the large window mode was used for the MOS cameras to avoid pile-up in Cyg OB2 #9 (which was then at its maximum). After pipeline processing, the data were filtered to keep only best-quality data (PATTERN of 0–12 for MOS and 0–4 for pn). Background flares were detected in several observations (Revs 0896, 0911, 1353, 1355, 2114, 3097). Only times with a count rate above 10 keV lower than 0.2–0.3 cts s⁻¹ (for MOS) or 0.4 cts s⁻¹ (for pn) were kept. A source detection was performed on each EPIC dataset using the task edetect_chain on the 0.3–2.0 keV (soft) and 2.0–10.0 keV (hard) energy bands and for a minimum log-likelihood of 10. This task searches for sources using a sliding box and determines the final source parameters from point spread function (PSF) fitting; the final count rates correspond to equivalent on-axis, full-PSF count rates (Table A.1). It may be noted in this table that the count rates of Schulte 12 in the off-axis observations are unreliable, as demonstrated by their much larger hardness ratios (due to uncorrected contamination by stray light).

We then extracted EPIC spectra of Schulte 12 using the task especget in circular regions of 35″ radius (to avoid nearby sources) centred on the Hipparcos positions (reported in Simbad). For ObsID 0793183001, the source appears so far off-axis that a larger ellipse was needed to extract it. On the contrary, in a few other cases the extraction radius had to be reduced for one of the EPIC cameras to avoid nearby CCD gaps. Whenever the core of the source fell in one such gap, the spectrum was not extracted to avoid large uncertainties. For the background, a circular region of 35″ radius was chosen in a region devoid of sources and as close as possible to the target. Dedicated ARF and RMF response files, which are used to calibrate the flux and energy axes, respectively, were also calculated by this task. EPIC spectra were grouped, with specgroup, to obtain an oversampling factor of five and to ensure that a minimum S/N of 3 (i.e. a minimum of 10 counts) was reached in each spectral bin of the background-corrected spectra.

2.2.2. Swift

At our request, the Cyg OB2 region was observed 18 times by the Neil Gehrels Swift Observatory. In addition, the same area was the subject of a short campaign in June–July 2015 and it serendipitously appears from time to time in the field of view of observations dedicated to PSR J2032+4127 (Petropoulou et al. 2018), 3EGJ2033+4118, or WR 144. All these XRT data taken in Photon Counting mode were retrieved from the HEASARC archive centre and were further processed locally using the XRT pipeline of HEASOFT v6.22.1 with calibrations v20170501. Because of the high extinction towards Cyg OB2, no optical loading is expected for XRT data of Schulte 12. The source spectra were extracted within Xselect using a circular region of 35″ radius; a nearby region devoid of any source was used for the background. The spectra were binned using grppha in a similar manner as the XMM-Newton spectra. The adequate RMF matrix from the calibration database was used whilst specific ARF response matrices were calculated for each dataset using xrtmkarf, considering the associated exposure map. Corrected count rates in the same energy bands as XMM-Newton were obtained for each observation from the UK online tool³, and they are reported in Table A.1.

UVOT photometry was taken in parallel with the X-ray data using a variety of filters. Following the recommendations of the Swift team⁴, magnitudes were extracted from Level2 files using an extraction region of 5″ radius for the source (and the same coordinates as in X-rays) while a nearby region of 15″ radius and devoid of sources was used as background. The resulting photometry is provided in the last column of Table A.1, where the reported errors correspond to the combination of systematic and statistical errors. No significant correlation is found between X-ray count rates and UV photometry.

2.2.3. Chandra

The massive stars of Cyg OB2 are X-ray bright, and therefore suffer from pile-up in usual Chandra-ACIS observations. However, Schulte 12 was also observed twice with Chandra-HETG (Table A.1). We reprocessed these data using CIAO 4.9 and CALDB 4.7.3. In the first dataset (ObsID = 2572), Schulte 12 appears far off-axis: its high-resolution spectrum is of low quality (Oskinova et al. 2017), but this position, leading to a larger PSF, allows its zero-order spectrum to avoid pile-up. This zero-order spectrum was therefore extracted using specextract in a circle of radius 20 px (9.8″) centred on the Simbad coordinates of the star and a nearby circle of 14.8″ radius was used as background region. The second exposure (ObsID = 16659) was pointed at Schulte 12, and therefore the high-resolution HEG and MEG spectra were extracted, combining the orders +1 and –1 for each grating, in addition to the zeroth-order spectrum. Specific detector response matrices (RMF and ARF) were calculated to obtain energy- and flux-calibrated spectra and a grouping similar to that used for XMM-Newton spectra was applied.

2.2.4. JAXA facilities

ASCA observed Cyg OB2 in April 1993 while Suzaku did so in December 2007. We downloaded these archival data from the HEASARC website and reduced them according to recommendations⁵.

The ASCA observation of Cyg OB2 consists of several datasets because of two different data modes (BRIGHT and BRIGHT2) and bit rates (high and medium). In all of them, Schulte 12 appears in the middle of stray light due to Cyg X-3. We extracted in Xselect the individual SIS-0 and SIS-1 spectra using a circle of 80″ radius, to avoid the emissions of the neighbouring stars, and considering the background from the surrounding annulus (with an external radius of 199″) to best compensate the stray light contamination. Data from high and medium bit rates were combined, but those from the two modes and two detectors were extracted separately, leading to four spectra in total. Following recommendations, RMF and ARF responses were calculated using the task SISRMG (for grades 0234) and the task ASCAARF (for a point source and simple case), respectively. A grouping similar to that used for XMM-Newton spectra was applied with GRPPHA. However, the final spectra were found to be extremely noisy and therefore unusable, so they were discarded from our analysis.

In Suzaku observations, we extracted in Xselect the individual XIS-0, XIS-1, and XIS-3 spectra of Schulte 12 in a circle of 80″ radius, to avoid the neighbouring emissions, and we used a nearby circle of 100″ radius for the background. Response files were automatically calculated in Xselect using the calibration database of June 2016. A grouping similar to that used for XMM-Newton spectra was applied with GRPPHA.

3. X-ray results

3.1. Light curves

As found previously with XMM-Newton (Cazorla et al. 2014), the Swift light curve in the total band clearly displays variations of a large amplitude, with rates appearing mostly between 0.04 and 0.10 cts s⁻¹ (Fig. 1). These changes do not appear correlated to hardness changes. To better characterise these variations in X-ray brightness, we applied a set of period-search algorithms: (1) the Fourier algorithm adapted to sparse/uneven datasets (Heck et al. 1985; Gosset et al. 2001, see also its recent rediscovery by Zechmeister & Kürster 2009 – these papers also note that the method of Scargle 1982, while popular, is not fully statistically correct), (2) three binned analyses of variances (Whittaker & Robinson 1944; Jurkevich 1971 which is identical, with no bin overlap, to the “pdm” method of Stellingwerf 1978, and Cuypers 1987 – which is identical to the “AOV” method of Schwarzenberg-Czerny 1989), (3) conditional entropy (Cincotta et al. 1999; Cincotta 1999, see also Graham et al. 2013), and (4) two different string length methods (Lafler & Kinman 1965; Renson 1978). The Fourier periodogram (top left panel of Fig. 2) clearly reveals a peak at P = 108.11 ± 0.32 d – the error here corresponds to the accuracy with which the period can be determined (i.e. a small fraction of the peak width). As expected from the spectral window (bottom part of the top left panel of Fig. 2), the peak is accompanied by secondary peaks, which are aliases of the main one and are not discussed further. The periodicity is confirmed by the other methods (that find P in the range 107.99–109.05 d), even though some of them, especially the string length methods, are less sensitive. If we divide the XMM-Newton EPIC-pn count rates by a factor of ten and make a period search on all XMM-Newton and Swift data combined, then a period of P = 107.99 ± 0.23 d is found by the Fourier method, with P = 107.64 − 108.45 d for the others (top right panel of Fig. 2).

Fig. 1. Swift light curve of Schulte 12.

Fig. 2. Top left: Fourier periodogram calculated from the Swift light curve, showing a clear peak at low frequencies. Significance levels of 0.1% and 0.01%, derived from tests (see text for details), are shown as dotted and dashed lines, respectively. Top right: comparison of period search results on the light curve composed of Swift data and scaled EPIC-pn data with different methods: the same peak (maxima for Fourier and AOV, minima for Renson, Lafler & Kinman, and conditional entropy) is detected in all cases and marked by a star. Bottom left: phase-folded X-ray light curve, using P = 107.99 d and T0 = 2 458 288.654. Our 2018 Swift monitoring campaign of Schulte 12 is shown in green and other high-quality (error on count rate < 0.008 cts s−1) Swift points in blue; EPIC-pn count rates divided by ten are shown as red stars; we note the good agreement between all data. Bottom right: same as previous, but with the lower-quality Swift data shown in dotted black.

The significance of this X-ray period can be readily estimated for the AOV method using the F statistics. The amplitude in the AOV periodogram corresponds to the ratio of the variance of the individual means in M phase bins and the sum over all phase bins of the variances about the mean within each bin. Amplitude therefore follows an F-distribution with (M − 1, N − M) degrees of freedom, with M the number of bins (chosen equal to 10) and N the number of data points (e.g. N = 96 for Swift data and 106 for Swift+XMM-Newton data). For the top-right AOV panel of Fig. 2, significance levels of < 0.01% are formally reached for amplitudes above 4.3: the observed peak amplitude therefore corresponds to a significance level well below 0.01%.

For the Fourier method, no analytical formula can be used to determine the significance level of the highest peak in case of uneven sampling, but the significance level can be estimated in other ways. First, we performed two sets of 10 000 Monte-Carlo simulations considering the observation times and drawing at random the individual count rates from a Gaussian distribution with the same mean for all simulations and a standard deviation equal to (1) the observational error on the count rate at each time or (2) the dispersion of the observed count rates. The first test probes the existence of the periodicity against a constant signal considering only the observational errors, while the second one tries to reproduce the variability range of the data assuming random variations, thereby testing a different null hypothesis. Secondly, we performed another test by shuffling the data 10 000 times: the count rate values were randomly associated with the observing times, again testing the periodicity considering the full range of the count rates. The 0.1% and 0.01% significance levels derived from these tests are shown on the top left panel of Fig. 2. With an amplitude of 0.014 cts s⁻¹, the observed peak is associated to a significance level below 0.01% for the first MC test and the shuffling test, and of 0.02% for the second MC simulation. With such low significance levels, the detection of the X-ray periodicity is thus undoubtly significant.

The phase-folded light curve (bottom panels of Fig. 2) further underlines the good agreement between Swift and XMM-Newton data. The confirmation of the period, in independent instruments and for data taken years apart, shows that the X-ray periodicity exists and is stable at least over the ∼5000 d interval covered by Swift and XMM-Newton.

3.2. Spectral fits

We fitted the spectra (Fig. A.1) under Xspec v12.9.1p, considering reference solar abundances from Asplund et al. (2009). As X-ray lines are clearly detected on both low- and high-resolution spectra (see also Oskinova et al. 2017), absorbed optically thin thermal-emission models of the type tbabs × phabs × ∑apec were considered, as is usual for massive stars. The first absorption component represents the interstellar contribution, fixed to 2.04 × 10²² cm⁻², a value calculated using the reddening of the target (E(B − V) = 3.33, Whittet 2015) and the conversion factor from Gudennavar et al. (2012).

We first fitted the XMM-Newton EPIC spectra, considering all available MOS and pn data simultaneously. Two thermal components were sufficient to achieve a good fit. Several combinations of absorption and temperature however yielded fits of similar quality (i.e. similar χ²), again a usual feature of the spectra of massive stars. We decided to fix the two temperatures to 0.75 and 1.9 keV, the most common values, lifting the absorption degeneracy, and repeated the fits. As we suspected, the fit quality was similar to what we had before. The same model was used for fitting Chandra and Suzaku spectra as well as the Swift spectra with the highest S/Ns (error on count rate below 0.08 cts s⁻¹, i.e. data shown as green and blue points in Fig. 2). The results of these fits are provided in Table 2 and Fig. 3. The fitted fluxes (and the normalization factors) closely follow the count rates, as expected (see Fig. 2). The hardness ratios do not vary significantly, but the absorption clearly increases when the observed flux is low.

Fig. 3. Results of the spectral fits: circumstellar absorption (top left), normalization factors (top right – the norm1 (resp. norm2) values appear in the upper (resp. lower) part, with filled (resp. empty) circles for XMM-Newton data), observed flux (bottom left), and ratio of unabsorbed hard and soft fluxes (bottom right). The results for XMM-Newton data are shown as black circles except for the stray light-contaminated cases displayed as red triangles; magenta stars indicate Chandra data; blue crosses the Suzaku observation; and Swift data are shown as empty green squares (shown for fluxes only). The phase folding is done as in Fig. 2.

One XMM-Newton point is clearly an outlier, with a stronger flux (Table 2) and a harder spectrum (Table A.1) than other XMM-Newton data. This point corresponds to the exposure taken in Rev. 3097 (ObsID 0793183001). As mentioned in Sect. 2, Schulte 12 appears far off-axis in this observation and within stray light from Cyg X-3. Despite a careful selection of background, it therefore appears that the stray light contamination cannot be totally eliminated. Therefore, we decided to discard these data from further consideration. We note that the spectra from the other off-axis observations (ObsID = 0801910201–601) appear less contaminated, with results usually within the errors (though they are far from perfect).

The Chandra and Swift data nicely agree with XMM-Newton, despite some remaining cross-calibration uncertainties (and a much larger noise for Swift), but the Suzaku results appear discrepant, especially in flux and norm₂ (the difference is less significant for the other parameters and the hardness ratio also seems in agreement). That Suzaku dataset presented neither an obvious stray light problem nor any reduction trouble. Furthermore, the closest XMM-Newton and Swift datasets, taken ∼200 d before and after the Suzaku dataset, respectively, agree well with the overall behaviour of Schulte 12; the origin of the Suzaku discrepancy therefore remains uncertain.

4. Optical results

4.1. Light-curve analysis

Gottlieb & Liller (1978) found Schulte 12 to be photometrically variable by up to 0.3 mag over decades while Laur et al. (2012) and Salas et al. (2015) reported changes on shorter timescales. Laur et al. (2012) found no clear periodicity, only a typical variation timescale of about 30 d. They also detected a potential trend in V − I colour, possibly linked to some spectral type variations though systematic errors could not be excluded. In contrast, Salas et al. (2015) detected significant changes with an amplitude in the I filter of 0.18 mag and with a period of 54.0 ± 0.1 d. This period is half the period found in X-rays, suggesting a common origin with two photometric events per X-ray period.

The picture is however more complex when one looks at larger datasets. Our recent monitoring of the star was performed at the same time as the spectroscopic campaign (Fig. 4). It clearly confirms the brightness variations of Schulte 12, with minima occurring at HJD ∼ 2 458 315, 2 458 370, and 2 458 460.

Fig. 4. Photometry recorded by one of the present authors (F.C., red): ASAS-SN (black), OMC (green), and NSVS (blue). The times of the Hermes and Carmenes spectroscopic observations are indicated by solid blue and dashed cyan vertical tick symbols, respectively.

Public photometric databases provide additional data. The noise affecting Hipparcos photometric measurements and Palomar Transient Finder (PTF) data renders them useless for detecting the small variations of Schulte 12, but significant variations are found in the Northern Survey for Variable Stars (NSVS)⁶ and in Integral-OMC⁷ observations, though their time coverage is scarce (Fig. 4).

ASAS-SN data⁸ (Shappee et al. 2014; Kochanek et al. 2017) in the g and V filters taken on similar dates to our monitoring agree well with it (as well as with OMC data). The ASAS-SN data further reveal previous minima around HJD ∼ 2 458 240 and possibly 2 457 950 (Fig. 4). While the separation between the two minima in Summer 2018 is about 54 d, the previous minima occurred ∼75 d and 365 d earlier than HJD  ∼  2 458 315. In addition, the minimum of Fall 2018 occurred 145 d later than that date, so the interval between minima does not seem to be constant. Using the X-ray ephemeris, these five minima correspond to phases of 0.86, 0.55, 0.24, 0.75, and 0.59 (considering the dates in chronological order). Only the third one (HJD  ∼  2 458 315) appears close to an X-ray minimum. Moreover, the photometric modulation seems to disappear at even earlier dates (HJD  ∼  2 457 000 − 800), with a long monotonic decrease detected between HJD  ∼  2 457 500 and HJD  ∼  2 457 700. No change in the ASAS-SN camera or observing strategy was performed however, meaning that this is not an instrumental effect. Salas et al. (2015) reported their detection of a periodic signal for observations taken in HJD ∼ 2 456 100 − 600. It therefore seems that the periodic signal may be transient.

Using the same methods as before, we performed a period search on the two largest datasets (ASAS-SN and OMC in V; see left panels of Fig. A.2). The periodogram based on ASAS-SN data only shows a peak at very low frequencies, at a position close to a peak in the spectral window: no significant periodicity is therefore detected (as could be expected from the discussion above). On the other hand, the OMC periodogram possesses a clear peak at about 87 d, but the scarce time sampling does not fully cover a single 87 d period (i.e. there are no set of 5–10 points taken within 87 d): it is therefore difficult to fully confirm it, but this signal could actually be related to the 54–108 d timescales. We also performed a period search on the most recent data only (last run of ASAS-SN data and our own photometric campaign; see right panels of Fig. A.2) since this is when the oscillations are most apparent. While they cover similar time intervals (see Fig. 4), though with different sampling, the two datasets do not yield the exact same periodograms, with the highest peak at 73.6 ± 2.2 d for ASAS-SN and at 48.4 ± 1.1 d or 130 ± 8 d for our own data. This confirms the probable absence of significant strict periodicity in the photometry, though significant variations on timescales of tens of days are clearly present.

4.2. Spectroscopic features

Apart from a large number of strong diffuse interstellar bands (DIBs, e.g. 5780, 5797, 6613 Å, Herbig 1995), the optical and near-IR spectrum of Schulte 12 contains prominent hydrogen and He I lines, but also a few metallic lines (Fig. 5). In particular, Paschen and Brackett hydrogen lines, He I λλ 5876, 6678, 7065 Å, Si II(+Mg II) λ 6347 Å, Si II λ 6371 Å, Ne I λ 6402 Å, N II λ 6482 Å, C II λλ 6578, 6583 Å, O I λλ 7772, 7774, 8446 Å, and N I λλ 8680–8718 Å appear in absorption, while He I λ 10 830 Å displays an incipient P Cygni profile, Hα displays a strong emission, and Fe II λ 7513, 9998 Å and Mg II λ 10952 Å show weaker emissions.

Fig. 5. Left: excerpts of the yellow-red and near-IR spectrum of Schulte 12, showing its main lines. Right: examples of line profile variations.

Criteria for the spectral classification of OB supergiants based on red spectra in the region around Hα and in the I-band were discussed by Negueruela et al. (2010). In the spectrum of Schulte 12, the rather narrow Paschen lines as well as the relative strength of the O I λ 8446 line compared to the neighbouring H I λ 8438 clearly indicate a mid-B supergiant star. On the other hand, the C II lines are less constraining since they are strong for giants and especially so for supergiants of spectral types between B1 and B5 (Walborn 1980; Negueruela et al. 2010). Although clearly present, the N II λ 6482 Å line is weak, suggesting a spectral type later than B2, which is confirmed by the presence of strong Si II lines (Walborn 1980). The presence of weak but definite N I lines in the region between 8680 and 8718 Å indicates a spectral type later than B3 but earlier than B6 (Negueruela et al. 2010). Comparing the equivalent widths (EWs) of Si II λ 6347 and Ne I λ 6402 (EW = 0.34 ± 0.02 Å and 0.16 ± 0.02 Å, respectively – see Table 3) with the results of Lennon et al. (1993, Figs. 21 and 24) further yields a spectral type around B5 with an uncertainty or a range of variation (see below) of about one subclass. The emission component of the He I λ 10 830 Å line is weaker than what is observed in very luminous (luminosity class Ia⁺) early B hypergiants (Groh et al. 2007). The weak Fe II λ 9998 Å emission line was seen in the spectra of HD 183143 (B7 Ia) and μ Sgr (B8 Iap) presented by Conti & Howarth (1999), but is absent in the spectrum of HD 225094 (B3 Ia) shown by the same authors. The above considerations suggest a spectral type B5 Ia for Schulte 12. This result is in reasonable agreement with the B3-4 Ia⁺ classification determined by Clark et al. (2012), except for the luminosity class. Whereas Clark et al. (2012) considered Schulte 12 to be a blue hypergiant, we favour a normal supergiant, which is also in line with the revised distance (see following section).

We compared the spectra taken close to minimum brightness (2018 May, July, and Sep) with those taken near maximum brightness (2018 June and mid-August). Several lines which are sensitive to spectral type (Lennon et al. 1993; Negueruela et al. 2010), and thus temperature, display changes of their strength that seem related to the brightness variations. For instance, the H I absorptions clearly appear stronger when the star is fainter and its V − I index is higher (Fig. A.3). A similar behaviour, though of smaller amplitude, is seen in Ne I λ 6402 Å, Si II λλ 6347, 6371 Å, and N I λ 8680 Å. Correlated changes are less obvious for the C II λλ 6578, 6583 Å, and N II λ 6482 Å lines, but they are less sensitive to temperature effects. The observed variations suggest changes of the spectral type between about B4 and B6. Such variations of the temperature are consistent with the V − I colour variations of about 0.05 mag that we found to be correlated with the brightness variations of the star.

Our dataset is composed of 31 high-resolution spectra (of which 19 also cover the near-IR wavelengths) obtained between 2010 and 2018, with a denser monitoring over the last year (one spectrum every ∼10 d). Comparing these spectra, it is clear that line profiles change from one observation to the next (Fig. 5). At some point, most lines seem to skew and even become double, a telltale sign of binarity. However, in the process, the depths of these lines do not change significantly, contrary to what is seen in spectroscopic binaries. Furthermore, other lines observed simultaneously display no sign of splitting, such as the Mg II and Fe II emissions which simply seem to shift. Finally, Hα and He I λ 10 830 Å show a more complex behaviour. The Hα profile is mostly composed of a broad emission with one or two peaks. In the latter case, both the relative intensity of the two peaks and their separation may vary. In addition, a shallow blueshifted absorption is sometimes (but not always) present. These variations are much more extreme than those reported by Chentsov et al. (2013) and Klochkova & Chentsov (2004). Finally, the He I λ 10 830 Å line profile is composed of three parts: a constant interstellar component (two deep and narrow absorption lines at velocities ∼0 and –30km s⁻¹), a broad emission component on the red wing, and a complex absorption on both wings. When overplotting all data, the broad emission appears constant, whereas the absorption features vary in number, depth, width, and position.

To further characterise the variability, we computed the temporal variance spectrum (TVS, Fullerton et al. 1996) for the main lines mentioned above (see a few examples in Fig. 6). To this aim, we combined the Hermes and Carmenes spectral sets, which can be done in view of their similar resolving powers and of the fact that the line widths are largely dominated by the star, not the instrument. The TVS confirms the presence of significant variations over a large range of velocities, larger than the −ν sin i  to  + ν sin i range (ν sin i = 38 km s⁻¹, see Clark et al. 2012). For all examined lines, the TVS displays a double-peaked structure. Such a profile is expected if the line centroid shifts back and forth in radial velocity. This can be the result of a genuine binary, pulsations, or co-rotating large-scale structures in the stellar wind. Another remarkable feature can be found when examining the symmetry of the TVS: the blue wing becomes progressively stronger when going from lines formed deep in the photosphere (such as the S II or C II lines which have roughly symmetric TVS profiles) to lines arising farther out (such as the hydrogen Paschen lines which have blue wings extending out to −150 km s⁻¹). This probably reflects increasing contributions from wind material to the absorption. The TVS of the Hα emission line reveals the largest and broadest variability (up to 15% of the continuum, between –200 km s⁻¹ and 150 km s⁻¹), with He I λ 10 830 Å being second (variations up to 9% in –150 km s⁻¹ and 100 km s⁻¹), but even the quieter Mg II emission exhibits changes by more than 1%.

Fig. 6. Mean spectrum and temporal variance spectrum for various lines in the spectrum of Schulte 12. The dotted lines in the TVS1/2 panels correspond to the 1% significance level of variability evaluated accounting for the S/N of the data. On the sides of the Mg II line, significant narrow variabilities are due to the imperfect correction of the telluric lines.

To characterise the lines, fitting profiles of a given shape (e.g. Gaussian) is a usual method. However, because of the variability of the profiles, it is quite difficult to obtain a good fit over the full profile – meaningful radial velocities can only be derived by fitting the bottom (resp. top) of the absorption (resp. emission) lines. Therefore, we decided to focus on line moments. Subtracting 1 from the normalized spectra, we calculated the zeroth-order moment M₀ = ∑(F_i − 1) which, multiplied by the wavelength step, yields the EW, the first-order moment M₁ = ∑(F_i − 1)×ν_i/∑(F_i − 1), which provides the centroid, and the second-order moment M₂ = ∑(F_i − 1)×(ν_i − M₁)²/∑(F_i − 1), whose square root gives the line width. Because of the noise, we avoided considering further orders. For each line, we limit the calculation to the wavelength interval where the line varies, using most of the profile while avoiding nearby lines. Table 3 displays the averaged values of these moments (⟨X⟩ = ∑X_i/N), along with their dispersions $(\sqrt{\sum {(X_i-⟨X⟩)}^2/(N-1)})$. All data were used in these averages, though a few noisy spectra (see Table 1) clearly add to the dispersion, especially for the bluer wavelengths. A very good agreement is found when comparing the two datasets, Carmenes and Hermes, demonstrating the reliability of the reduction. The diffuse interstellar bands (DIBs) enable us to check the (small) “natural” noise around the moment values, due to the number of collected photons and the quality of the reduction (wavelength calibration, normalization, etc.). This provides a threshold against which variability can be detected. Equivalent widths and widths generally are rather stable except for the two lines strongly affected by wind, Hα and He I λ 10 830 Å. It is interesting to note that the dispersions in EW and width reach their smallest values not only for DIBs but also for Mg II and Fe II lines: this confirms that these stellar lines vary very little in shape. The centroids are however found to vary for all stellar lines, with the largest dispersions found for the wind lines and the smallest (again) for DIBs, and then for Mg II and Fe II emission lines.

To clarify the behaviour of the lines, we compared the individual centroids of different lines in Fig. 7. Except for the wind lines, strong correlations between centroids are found. But the devil lies in the details. Indeed, the amplitudes of variations fall into two groups: the centroids of He I λλ 5876, 6678, 7065 Å, Si II λλ 6347, 6371 Å, C II λλ 6578, 6583 Å, and H I λλ 8598, 8863, 10 049, 15 556 Å (all absorption lines) vary by about 25 km s⁻¹ peak-to-peak, whereas those of Fe II λλ 7513, 9998 Å, and Mg II λ 10 952 Å emissions only change between –20 and –30 km s⁻¹. Indeed, a unity slope can be used for all cases except when comparing centroids from the first group with those of the second (see bottom right panel of Fig. 7). However, it is also clear that even the variations of Fe II and Mg II are much larger than that observed for DIBs.

Fig. 7. Comparisons between the centroid of He I λ 6678 Å and the centroids of other lines detected in the spectrum of Schulte 12 (He I lines on top left, H I lines on top right, and metallic absorption and emission lines on bottom left and right, respectively). Filled symbols correspond to Hermes data, open symbols and crosses to Carmenes data; the diagonal yields the one-to-one correlation. For comparison, the case of a DIB is added to the last panel.

Rivinius et al. (1997) reported the results of an intense monitoring of three early-B hypergiants. All three stars displayed significant velocity variations of their photospheric lines. For instance, in the case of ζ ¹ Sco, Rivinius et al. (1997) observed a pulsation-like pattern with a peak-to-peak amplitude of 20 km s⁻¹ on timescales near 10–15 d that gave way after some time to a more random variation. Lines originating from outer layers of the photosphere were found to display a delayed variation pattern compared to lines formed in deeper layers. In addition, we may note that photometric variations with an amplitude near 0.1 mag were observed for ζ ¹ Sco on the same timescales (Burki et al. 1982). A similar situation is observed in the case of Schulte 12 (though with a longer variation timescale, 50–100 d): the velocities of He I λ 6678 Å and H I λ 8863 Å lines display the same qualitative behaviour as those of C II λλ 6578, 6583 Å lines, but with a slight delay of the order of 3–5 d (see Fig. 8).

Fig. 8. Evolution of the photometry of Schulte 12 and of the centroids of some prominent lines in its spectrum during the Carmenes campaign in Summer 2018.

Finally, we performed period searches on the line profiles using the modified Fourier algorithm (Heck et al. 1985) and the same intervals as for the moment calculation. To this aim, we interpolated all profiles to a common velocity grid with steps of 0.2 km s⁻¹, and then calculated periodograms at each velocity considering the same velocity intervals as used for moment calculations. These individual periodograms were then averaged to get a single periodogram per line (Fig. 9). While no frequency stands out, multiple peaks often appear near 0.015–0.025 d⁻¹ for absorption lines, hinting at long-term variability on timescales of 40–65 days, that is, close to about half the X-ray period. Indeed, spectra taken only a few days apart (2017 May 24 and 29) show no significant difference while variations are detected on longer timescales, typically tens of days (see Figs. 4 and 8 – e.g. in the latter figure, the times of the two minimum velocities are separated by ∼50 d).

Fig. 9. Periodograms calculated on the line profiles (see text) for H I λλ 8598, 8663, 10 049, 15 556 Å, He I λλ 5876, 6678, 7065 Å, C II λλ 6578, 6583 Å, O I λλ 7772, 7774 Å, Mg II λ 10 952 Å, and the DIB at 6613 Å for comparison. The vertical dotted line shows the frequency corresponding to the 108 d period.

5. Discussion

The previous sections led to several results. Firstly, at high energies, the flux varies by a factor of about two and a period of 108 d is detected for these changes. The minimum flux corresponds to a time of lower intrinsic emission and increased absorption. Secondly, we confirmed the presence of photometric and spectroscopic changes in the star, which we further constrained. However, no regular clock seems present in the optical/near-IR domain, a simple timescale of 50–100 d being found. We now attempt to interpret these contrasting features.

5.1. The evolutionary status of Schulte 12

The second data release of Gaia (Gaia Collaboration 2018) allowed for the first time a direct measurement of the parallax of Schulte 12 to be made. This led to a surprising value of 1.18 ±0.13 mas. Indeed, this must be compared to the DR2 parallaxes of the other three main objects of the Cyg OB2 association, which are compatible with the usual distance used for Cyg OB2 (e.g. d∼1.75 kpc in Clark et al. 2012): 0.64 ± 0.06 mas for Cyg OB2 #5, 0.62 ± 0.04 mas for Cyg OB2 #8A, and 0.60 ± 0.03 mas for Cyg OB2 #9. Schulte 12 thus appears twice closer than the stars with which it is often associated. Proper motions also appear quite different, with values of −1.9 ± 0.2 mas in RA and −3.3 ± 0.2 mas in DEC for Schulte 12 and values of minus 2.7–3.1 mas in RA and minus 4.1–4.7 mas in DEC for the others. This actually confirms early doubts on the membership of Schulte 12 to Cyg OB2 (Walborn 1973).

One may however wonder whether or not the presence of visual companions to Schulte 12 could impair the Gaia measurements. In this context, it is important to recall that all four X-ray bright massive stars of Cyg OB2 are actually multiple, with very diverse properties: Cyg OB2 #8A and 9 consist of spectroscopic binaries with periods of 22 d and 2.4 yr, respectively, while Cyg OB2 #5 is thought to have four components with periods ranging from 6.6 d to 9000 yr. If companions were modifying the parallax and proper motions of Schulte 12, it would be very surprising to find that the data of the other three systems agree – but they do. Therefore the parallax difference of Schulte 12 is very probably real. In this context, one may also argue that the angular diameter of Schulte 12 amounts to 1.31 mas, i.e. comparable to the parallax, as is the case of other blue supergiants (ζ¹ Sco, HD 190603, BP Cru, HD 80077, HD 168607, and to a lesser extent P Cygni and HD 168625). Shifts in the photocentre of the star (e.g. due to non-radial pulsations or spots on the stellar surface; see below) could thus affect the parallax. However, we note that the absence of a strong, strictly periodic photometric modulation does not support this scenario. Furthermore, in view of the rather small relative error on the Gaia-DR2 parallax, it would require an extraordinary coincidence for the shifts of the photocentre to exactly mimic the effect of an annual parallax in the Gaia data. We thus consider the revision of the distance as robust.

This new distance (${840}_{755}^{945}$ pc in Bailer-Jones et al. 2018) drastically changes the accepted physical properties of Schulte 12: scaling the values of Clark et al. (2012) yields a luminosity of log[L_BOL/L_⊙] = 5.64, a radius of 118 R_⊙, and a spectroscopic mass of 25 M_⊙. The mass-loss rate proposed by Clark et al. (2012) or Morford et al. (2016) would also be reduced by about a factor of four. Schulte 12 now appears as a normal supergiant, without any luminosity problem (a difficulty noted by Clark et al. 2012).

It is interesting to compare Schulte 12 with the LBV R 71 in the Large Magellanic Cloud which is at the lower luminosity end of the classical LBV strip (Mehner et al. 2017). Both stars share similar properties (Mehner et al. 2017; Clark et al. 2012, and this work): they have similar spectroscopic masses (27 M_⊙ for R 71 and 26 M_⊙ for Schulte 12), temperatures (15 500 ± 500 K for R 71 in quiescence and $13{700}_{-500}^{+800}$ K for Schulte 12), radii (107 R_⊙ for R 71 in quiescence and 118 R_⊙ for Schulte 12), and luminosities (6.0 10⁵ L_⊙ for R 71 in quiescence and 4.4 10⁵ L_⊙ for Schulte 12). The slight luminosity difference is sufficient to keep Schulte 12 outside of the LBV strip, but actually close to it, while R 71 lies inside of it. Evidence that Schulte 12 is less evolved than R 71 comes from the chemical composition of its atmosphere. Mehner et al. (2017) reported a strong enrichment of the atmosphere of R 71 with CNO processed material (with number ratios of He/H = 0.20, O/N = 0.08 and C/N = 0.05). For Schulte 12, Clark et al. (2012) derived a significantly weaker enrichment (He/H = 0.1, O/N = 1.48, C/N = 0.31).

5.2. Binarity

As already mentioned in Sect. 4.2, the changing profile of the absorption lines cannot be due to line blending as the profile broadens, skews, or becomes narrower without much changes in line depth (Fig. 5). However, the Mg II and Fe II emissions undergo simple velocity shifts, which could be interpreted in this framework. Combining the centroid velocities of the Mg II and Fe II emission lines, we computed using LOSP⁹ a highly preliminary orbital solution. In this calculation, we assume that the 108 d X-ray period is linked to orbital motion. Freeing the period does not significantly change the quality of the fit or the values of the orbital parameters. The results are given in Table 4 and illustrated in Fig. 10. The periastron passage of this tentative orbit occurs ∼10 d (or Δ(ϕ) = 0.1) after the X-ray maximum and Δ(ϕ) = 0.25 before the X-ray minimum - in colliding-wind binaries however the periastron usually corresponds quite closely to an X-ray extremum.

Fig. 10. Tentative orbital solution of Schulte 12 for a fixed orbital period of 107.99 d (middle column of Table 4). The average velocities of Mg II and Fe II emissions in the Carmenes dataset are shown by blue dots and the Hermes velocities of Fe II λ 7513 Å (not used to obtain the orbital solution) by open red squares. We note that the phase used here corresponds to the binary ephemeris (Table 4), not to the X-ray ephemeris used for Figs. 2 and 3.

If we assume Schulte 12 to be a binary system, we need to ask ourselves whether or not the spectroscopic radius of the B supergiant is consistent with the radius of its Roche lobe. For this purpose, we use the formula of Eggleton (1983) to express the dependence of the Roche lobe radius on the mass-ratio $q=\frac{m_1}{m_2}$ as follows:

$$\begin{array}{c}\hfill \frac{R_{\mathrm{RL}}}{R_{\odot }}=\frac{2.06q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})}{\left(\frac{1+q}{q}\right)}^{1/3}{\left(\frac{m_1}{M_{\odot }}\right)}^{1/3}{\left(\frac{P}{\mathrm{day}}\right)}^{2/3}·\end{array}$$(1)

If we adopt a period of 108 d and a spectroscopic mass of m₁ = 25 M_⊙ for the B supergiant, we find that the radius of the Roche lobe of the B supergiant exceeds 127 R_⊙ whatever the value of q. Therefore, if the orbit were circular, the spectroscopic radius of 118 R_⊙ would fit entirely into the Roche lobe. The B supergiant would fill less than 80% of the volume of its Roche lobe, regardless of the mass of its companion. However, this is valid for a circular orbit, not in the presence of a significant orbital eccentricity as found here. Assuming that in a long-period eccentric binary the size of the Roche lobe varies in proportion to the instantaneous orbital separation, we find that for an eccentricity of about 0.5 (Table 4) the B supergiant would significantly overflow its Roche lobe at phases near periastron, regardless of the value of the mass of the companion.

Nevertheless, there is a significant dispersion of the data points around the best-fit curve, and moreover, plotting the additional Fe II λ 7513 Å velocities measured on Hermes spectra does not reveal a coherent behaviour (Fig. 10). It must be underlined here once again that Hermes and Carmenes spectra agree well: the spectra taken in May 2017 with both instruments at only a 5 d time difference appear very similar while the DIB centroids coincide within errors (see Table 3). Therefore, no problem in the instrument or during the analysis can be blamed for the difference observed in Fig. 10. In fact, no orbital solution using the full Fe II λ 7513 Å velocity dataset can be derived as these data simply cannot be reconciled with a 108 d period. Combined to the partial correlation of the emission line velocities with those of absorptions (Figs. 7 and 8), this suggests that we are not witnessing pure orbital motion. If Schulte 12 is a binary, we have not yet found evidence for it in the RVs.

Finally, we examined the spectrum of Schulte 12 in search of lines moving in the opposite direction compared to Mg II or of lines with a different ionization (e.g. He II), as they could then be attributed to a companion. None were found however. With the revised distance of Schulte 12, its absolute V-band magnitude amounts to −8.26. We estimate that our spectroscopic data should enable us to detect the spectroscopic signature of an O-type companion down to 2.5 mag fainter than the B supergiant. Considering the typical parameters of O-type stars from the calibration of Martins et al. (2005), we conclude that a main-sequence O-type star would remain entirely hidden, regardless of its spectral type. For O-type giants, only objects with a spectral type earlier than O5 would show up in the combined spectrum. Finally, O supergiants should be detectable in a combined spectrum regardless of the spectral types. The non-detection of even a massive companion would therefore not be entirely surprising.

5.3. Pulsational activity

The erratic variations of the centroids and profiles of the absorption lines suggest that the velocities of these lines are dominated by effects other than genuine orbital motion. One explanation for these variations could be non-radial pulsations of the B supergiant; but is Schulte 12 expected to be pulsationally unstable?

With the revised physical properties, the position of Schulte 12 in the Hertzsprung–Russell diagram (Fig. 11) drastically changes. The star now appears relatively far away from the empirical Humphreys–Davidson limit (Humphreys & Davidson 1994) and thus from the region where strange-mode oscillations are expected (Kiriakidis et al. 1993). Maeder (1980) suggested that intermediate term variability of B–G supergiants could stem from non-radial gravity-mode pulsations connected to a large outer convection zone. Indeed, Schulte 12 is located within the part of the Hertzsprung–Russell diagram where oscillatory convection modes with ℓ ≤ 2 are expected (Saio 2011). If we attribute the photometric variability of Schulte 12 to non-radial pulsations, they must indeed be of low ℓ to explain the relatively large observed amplitude of the variations. For a 30 M_⊙ star, Saio (2011) predicts periods of excited ℓ = 1 and ℓ = 2 modes of about 70 and 60 days, respectively. These values are reasonably close to the timescale of the variations of the absorption lines and the photometry of Schulte 12. Furthermore, the pulsation modes are associated to subphotospheric convection which is known to be unstable, hence the pulsational activity could be intermittent, and they have growth timescales comparable to their period, suggesting that such (dis)appearance could occur quite quickly, in line with our observations of Schulte 12.

Fig. 11. Revised position of Schulte 12 in the Hertzsprung-Russell diagram adopting the distance corresponding to the Gaia-DR2 parallax. The evolutionary tracks are taken from Ekström et al. (2012) for solar metallicity and no rotation. The continuous red line corresponds to the empirical Humphreys–Davidson limit (Humphreys & Davidson 1994). Various instability boundaries are shown according to Saio (2011). The dotted red curve indicates the instability boundary of low-order radial and non-radial p-mode oscillations. The short-dashed magenta contour indicates the region where non-radial g-modes with ℓ ≤ 2 are excited by the iron opacity bump. Finally, the cyan long-dashed contour indicates the region where oscillatory convection modes with ℓ ≤ 2 are excited (Saio 2011).

The line profile variability (LPV) of absorption lines due to non-radial pulsations depends on the ratio k between the azimuthal and radial pulsation velocities (Kambe & Osaki 1988). The fact that the strongest LPV signatures are observed in the cores of the lines (see Fig. 6), suggests that k ≪ 1. Schulte 12 lies also at the edge of the zone where low-order (n = 0 or n = 1) radial pressure modes are expected to be excited according to Jeffery & Saio (2016). However, the corresponding fundamental (n = 0) radial modes should have periods of order 10⁶ s, that is about 12 days, which is much shorter than the observed periods. It seems therefore unlikely that the observed variability of Schulte 12 could arise from such radial modes.

At this stage, it is once more interesting to compare Schulte 12 with the LBV R 71. During quiescence, R 71 displays microvariations of 0.1 mag amplitude on timescales of 14–100 days (Mehner et al. 2017, and references therein). This is qualitatively similar to what we observe in Schulte 12, probably indicating that both stars are subject to similar pulsation mechanisms. During the currently ongoing outburst of R 71, a new kind of nearly periodic photometric variation appeared, consisting of oscillations of 0.2 mag amplitude on a timescale of about 445 days (Walborn et al. 2017; Mehner et al. 2017). The star R 71 currently displays a late F supergiant spectrum (Walborn et al. 2017; Mehner et al. 2017), very different from that of Schulte 12. However, an interesting point concerns the spectroscopic variability that occurs on a timescale twice that of the above-mentioned 445 day photometric cycle. Indeed, the absorption lines vary between single and double morphologies (Mehner et al. 2017; Walborn et al. 2017), in a way similar to what we observe for Schulte 12. By analogy with RV Tau variables, Mehner et al. (2017) attribute the line doubling to shock waves triggered by the pulsations and propagating through the atmosphere of R 71. Whether a similar explanation holds in our case is unclear, as the atmospheric structure of the B5 Ia supergiant Schulte 12 is rather different from that of a late-F supergiant.

5.4. The X-ray clock

The X-ray data for Schulte 12 display a clear periodicity of 108 d. This is incompatible with the putative periods of the visual companions (see Introduction), and therefore they do not play a role in this phenomenon. The question then arises regarding the origin of the high-energy emission.

The folded X-ray light curve of Schulte 12 is strongly reminiscent of that of the massive binary WR 21a, where X-rays arise in the collision between the stellar winds of the two stars (Gosset & Nazé 2016). As in Schulte 12, the count rate remains relatively constant most of the time (or slightly increases), then a stronger increase occurs, followed by a sharp drop and a longer recovery. In WR 21a, the increase was associated to the typical 1/D increase in flux found in adiabatic wind–wind collisions as stars approach periastron (Stevens et al. 1992). The subsequent drop combined two effects: an initial increase of absorption as the companion plunged towards the densest wind layers at periastron, then a shock collapse followed by a slow recovery (Gosset & Nazé 2016). The two X-ray light curves thus appear very similar, except for the amplitude of the X-ray variation (a factor of five between minimum and maximum in WR21a rather than two in Schulte 12) but this can easily be accounted for by different physical and orbital parameters, as the X-ray properties are very sensitive to them. Further supporting this scenario is the presence of very hot plasma (demonstrated by the need of kT₂ = 1.9 keV and the presence of the iron line at 6.7 keV), as well as the strong X-ray luminosity of Schulte 12 (the flux in the 0.5–10 keV energy band, corrected for interstellar absorption, lies between 5.6 and 11.9 × 10⁻¹² erg cm⁻² s⁻¹, corresponding to 4.7 − 10 × 10³² erg s⁻¹ or log[L_X/L_BOL] = − 6.6 to –6.2). Both are common features in wind–wind collisions.

However, to produce a significant wind–wind collision, the companion of Schulte 12 ought to be a massive star of spectral-type O, early B, or Wolf–Rayet with a fast wind and in an eccentric orbit (to produce the 1/D effect). From the spectral analysis in the optical, such a companion cannot be an O-type supergiant or early giant, but other possibilities remain open in view of the brightness of Schulte 12. Nevertheless, as the spectroscopic mass of Schulte 12 is rather modest, the required massive OBWR companion would then lead to significant orbital motion for Schulte 12, except if the orbital inclination is low. In that case however an increase of the absorbing column towards a putative wind–wind collision zone could only occur at periastron, contrary to observations. In fact, as mentioned above, we certainly cannot claim to have detected the binarity nature of Schulte 12, hence the wind–wind collision scenario appears quite unlikely.

As the optical analysis showed, non-radial pulsations are most probably at work in Schulte 12. Could the X-ray light curve then be driven by pulsation? X-ray pulsations occurring in phase with the optical pulsations have indeed been observed for a few β Cep-type pulsators (Oskinova et al. 2014; Nazé 2015; Cazorla et al. 2014). There are however a number of caveats with this hypothesis. Firstly, the X-ray light curves of ξ¹ CMa (Oskinova et al. 2014; Nazé 2015) and HD 44743 (Cazorla et al. 2014) are nearly sine waves, whereas the X-ray light curve of Schulte 12 is far more asymmetric. It must also be mentioned that the (X-ray) pulsation periods of the two β Cep stars are much shorter (4.9 and 6.0 h for ξ¹ CMa and HD 44743, respectively) than the period of the X-ray light curve of Schulte 12. Secondly, the peak-to-peak amplitudes of the modulation of the X-ray fluxes in ξ ¹ CMa and HD 44743 are less than 10% of the mean flux. In the case of Schulte 12 this relative peak-to-peak amplitude reaches 70% of the mean flux. Thirdly, the modulation of the X-ray fluxes of the two β Cep stars occurs exactly on the same period as the modulation of their optical light curve. In the case of Schulte 12, the X-ray flux is modulated on a period roughly twice that of the optical variations (when they exist). Lastly, with log[L_X/L_bol] = − 6.7 and −7.6 for ξ ¹ CMa and HD 44743, respectively, the level of X-ray emission of the two β Cep stars is very close to (or slightly below) the typical one of massive stars, while Schulte 12 clearly appears overluminous in X-rays. The two β Cep stars also display the typical soft X-ray emission (kT <  1 keV, Nazé et al. 2014; Cazorla et al. 2014). There is therefore no observational evidence in favour of this scenario, and a detailed modelling for the peculiar pulsations of Schulte 12 would be required to demonstrate their influence on the X-ray emission.

Another candidate for a regular clock is rotation, as the 108 d period is compatible with the rotation estimates of Schulte 12. Clark et al. (2012) derived ν sin i = 38 km s⁻¹ for Schulte 12: adopting a radius of 118 R_⊙ leads to an upper limit on the rotational period of 158 days. Likewise, adopting the spectroscopic mass, radius, and luminosity, we infer a lower limit on the rotational period of 41 days, assuming the equatorial rotation velocity reaches the break-up velocity. Any observed periodicity in this rather wide range could thus possibly stem from rotational modulation. Two possible scenarios can be envisaged in this context: spots on the stellar surface leading to co-rotating interaction regions (CIRs, Cranmer & Owocki 1996; Lobel & Blomme 2008), and a magnetically confined stellar wind (Babel & Montmerle 1997; ud-Doula & Nazé 2016).

Considering the first case, it should be noted that variability of the Hα emission line is ubiquitous among OBA-type supergiants such as Schulte 12 (e.g. Fullerton et al. 1997; Kaufer 1999; Morel et al. 2004). Based on the similarity between the observed recurrence timescales and the estimated rotation periods (or an integer fraction of them), this phenomenon is frequently attributed to co-rotating large-scale structures in the wind (Fullerton et al. 1997). Investigating a sample of 22 OB supergiants with spectral types in the range O7.5 to B5, where Schulte 12 lies, Morel et al. (2004) found significant Hα variability for all of them. These latter authors further reported the existence of a cyclical pattern of variability in the Hα profile of nine stars of their sample. They also noted the existence of periodicities ranging between 1 and 27 days in the photometric data of nine stars. Hence the presence of CIRs in Schulte 12 would not be a surprising occurrence. In this context, several cases of rotational modulation of the X-ray flux of O-stars have been reported in recent years and were attributed to such large-scale structures (Oskinova et al. 2001; Nazé et al. 2013, 2018b; Massa et al. 2014; Rauw et al. 2015b). However, their variation timescales are shorter (1.8–4.2 d) than for Schulte 12 and the peak-to-peak amplitude of their modulations is ∼10% of the mean flux, which is much smaller than for Schulte 12. Furthermore, the X-ray brightness of these stars are in agreement with expectations for massive stars and their plasma temperature remains low, kT <  1 keV (Nazé 2009; Rauw et al. 2015b; Nazé & Motch 2018), contrasting with Schulte 12. Finally, while the optical light curve of Schulte 12 is strongly epoch-dependent, the photometric variations of such stars remain relatively coherent over timescales of years (see Morel et al. 2004 for a sample of B-supergiants and Ramiaramanantsoa et al. 2014, 2018 for ξ Per and ζ Pup, respectively). There is thus poor observational support for this scenario.

In the second case, X-ray overluminosities, large X-ray modulations, and large plasma temperatures may occur (Nazé et al. 2014) but this observational picture was mostly drawn from magnetic main sequence stars and therefore its applicability to supergiants remains to be demonstrated. Magnetic fields are seldom detected in massive supergiants, not only because of a (expected) low incidence rate but also for technical reasons. Indeed, as the stars evolve, their radii increase and their magnetic fields, if descending from the main-sequence fields, should be much smaller and therefore more difficult to detect. For these reasons, only a few magnetic detections of supergiants have been achieved until now (Blazère et al. 2015; Martin et al. 2018) and for ζ Ori A (O9.5I) the X-ray emission appears in line with the “typical” properties of massive stars, that is it seems unaffected by the presence of the magnetic field (Nazé et al. 2014). The applicability of this scenario is therefore suspended to a detailed investigation of the wind properties of Schulte 12 combined with the (positive) results of a spectropolarimetric campaign. In addition, our optical spectroscopy did not reveal a strict periodic modulation of the Hα emission line contrary to observations and expectations of magnetically confined winds (Townsend et al. 2005; Sundqvist et al. 2012; ud-Doula et al. 2013).

A last, tentative possibility would be sporadic accretion onto a companion of the B-star in an eccentric orbit. The pulsations could then be unrelated or excited through tidal interaction at periastron. However, the X-ray emission of Schulte 12 is thermal in nature, thus different from the typical non-thermal emission associated to accretion (though some peculiar accretion regimes may produce thermal emission). Furthermore, the X-ray variations concern not only the intrinsic strength of the emission but also the local absorption and therefore they are not in line with a simple additional emission appearing and disappearing at some phases.

6. Conclusions

We performed a multiwavelength campaign dedicated to Schulte 12, which we complemented with archival data. This unique dataset sheds new light on the properties of this star.

Firstly, an X-ray clock is present in the system: a period of 108 d, with a minimum flux occurring ∼35 d later than the maximum flux. This minimum corresponds to an increase in absorption along the line of sight and a decrease in the intrinsic emission of Schulte 12. The same phased light curve shape is found (1) from data scattered over a 14 year timescale (or 47 cycles, demonstrating the stability of the signal) as well as in a dedicated monitoring covering a full cycle (enabling us to reject aliases or incorrect period identifications); and (2) from fully independent instruments onboard the Swift and XMM-Newton telescopes, which observed Schulte 12 with very different temporal sampling. The period is detected by different period-search methods dedicated to uneven samplings; its significance level (modified Fourier, AOV) is ≤0.02%. Therefore the periodicity does not appear to be spurious.

Secondly, we confirmed the spectral type of Schulte 12 as well as previous reports of its photometric and spectroscopic variability, which we further characterised. While not strictly periodic, their timescale amounts to tens of days, in agreement with the X-ray period (or rather half of it). In the visible spectrum, all absorption lines of hydrogen, helium, and metals behave in a similar way. Their strengths (EWs) and widths do not vary much, but their centroids and shapes suffer from significant changes – for example velocities cover a range of about 25 km s⁻¹, while lines alternate between blueward-skewed, rightward-skewed, symmetric, and doubling profiles. The velocities of H and He I lines appear to lag by a few days behind those of metals, probably due to different formation depths. Emission lines of Mg II and Fe II appear less variable, with profiles simply shifting by 10 km s⁻¹, but these variations, which are larger than the noise level, do not seem to indicate binary motion. Finally, Hα and He I λ 10 830 Å, the two lines most affected by wind emission, display spectacular changes, of an amplitude larger than reported before. The behaviour of these lines does not seem obviously linked to that of the other lines and they do not exhibit strict periodicities either.

The Gaia data reveal that Schulte 12 is actually twice closer than other massive stars in Cyg OB2. This led to a severe revision of its stellar parameters: Schulte 12 is no longer a hypergiant, but a normal supergiant, and is no longer close to the empirical Humphreys–Davidson instability limit; it is however located in an instability zone linked to convection modes for which periods of tens of days are expected. This certainly explains the changes in the absorption line profiles, though the absence of strict periodicity implies intermittent activation of the oscillations.

Regular modulations of the X-ray emission of massive stars have been observed in several cases, but none seem to strictly apply to Schulte 12. The amplitude and shape of the X-ray light curve, the presence of a hot plasma, and the significant X-ray overluminosity recall the characteristics of colliding-wind binaries. However, the lack of any manifest binary signature in the optical/near-IR spectrum seems to exclude this possibility. In a few cases, pulsations have been linked to X-ray activity, but without overluminosity or high-temperature plasma and with shorter timescales and smaller variability amplitude than observed for Schulte 12. Rotational modulations of the X-ray light curves also exist. When attributed to CIRs, they again lack the overluminosity, hot plasma, and large variability amplitude observed in Schulte 12. When linked to the presence of magnetic field, such properties are encountered, but only for main sequence stars. The origin of the X-ray clock therefore remains mysterious. Further investigation, for example of UV wind lines, is now required.

Acknowledgments

Y.N. and G.R. acknowledge support from the Fonds National de la Recherche Scientifique (Belgium), the Communauté Française de Belgique, the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme (contract XMaS). S.C. acknowledges support through DFG projects SCH 1382/2-1 and SCHM 1032/66-1. We thank the Swift team for their kind assistance. This project has received funding from the European Union’s Horizon 2020 research and innovation programme via the OPTICON project under grant agreement No 730890. This material reflects only the authors views and the Commission is not liable for any use that may be made of the information contained therein. We are grateful to the Calar Alto service mode observers (M. Moreno, B. Arroyo, A. Guijarro, I. Hermelo, R.P. Hedrosa, P. Martin, G. Bergond & A. Fernández) for their help in collecting the Carmenes data. Some results were obtained with the Mercator telescope, operated by the Flemish Community on the island of La Palma at the Spain Observatory del Roche de los Muchachos of the Instituto de Astrofisica de Canarias. Mercator and Hermes are supported by the Funds for Scientific Research of Flanders (FWO), the Research Council of KU Leuven, the Fonds National de Recherche Scientifique (FNRS), the Royal Observatory of Belgium, the Observatoire de Genève, and the Thüringer Landessternwarte Tautenburg. L.M. thanks all the observers who collected the data with Hermes. ADS and CDS were used for preparing this document.

References

Appendix A: Complementary table and images

Fig. A.1. X-ray spectrum of Schulte 12 (left: Swift data of observation 00032767003, right: XMM-Newton data taken on rev. 3176 with pn in green, MOS1 in black, and MOS2 in red) along with their best-fit models (see Table 2) and the residuals.

Fig. A.2. Results from period searches on the visible light curves. From left to right: results for OMC data, ASAS-SN data (all runs and last run only), and our own photometry. In each panel, from top to bottom: periodogram resulting from the Fourier method (with its spectral window), the AOV method, and the conditional entropy method (see text for details).

Fig. A.3. Variation of the line strengths and the photometry in 2018. Formal errors on EWs are similar to the symbol size. The dashed cyan and dotted black vertical lines indicate times of maximum and minimum brightness, respectively.

All Tables

All Figures

	Fig. 1. Swift light curve of Schulte 12.
In the text

Fig. 2. Top left: Fourier periodogram calculated from the Swift light curve, showing a clear peak at low frequencies. Significance levels of 0.1% and 0.01%, derived from tests (see text for details), are shown as dotted and dashed lines, respectively. Top right: comparison of period search results on the light curve composed of Swift data and scaled EPIC-pn data with different methods: the same peak (maxima for Fourier and AOV, minima for Renson, Lafler & Kinman, and conditional entropy) is detected in all cases and marked by a star. Bottom left: phase-folded X-ray light curve, using P = 107.99 d and T0 = 2 458 288.654. Our 2018 Swift monitoring campaign of Schulte 12 is shown in green and other high-quality (error on count rate < 0.008 cts s−1) Swift points in blue; EPIC-pn count rates divided by ten are shown as red stars; we note the good agreement between all data. Bottom right: same as previous, but with the lower-quality Swift data shown in dotted black.

In the text

Fig. 3. Results of the spectral fits: circumstellar absorption (top left), normalization factors (top right – the norm1 (resp. norm2) values appear in the upper (resp. lower) part, with filled (resp. empty) circles for XMM-Newton data), observed flux (bottom left), and ratio of unabsorbed hard and soft fluxes (bottom right). The results for XMM-Newton data are shown as black circles except for the stray light-contaminated cases displayed as red triangles; magenta stars indicate Chandra data; blue crosses the Suzaku observation; and Swift data are shown as empty green squares (shown for fluxes only). The phase folding is done as in Fig. 2.

In the text

	Fig. 4. Photometry recorded by one of the present authors (F.C., red): ASAS-SN (black), OMC (green), and NSVS (blue). The times of the Hermes and Carmenes spectroscopic observations are indicated by solid blue and dashed cyan vertical tick symbols, respectively.
In the text

	Fig. 5. Left: excerpts of the yellow-red and near-IR spectrum of Schulte 12, showing its main lines. Right: examples of line profile variations.
In the text

	Fig. 6. Mean spectrum and temporal variance spectrum for various lines in the spectrum of Schulte 12. The dotted lines in the TVS1/2 panels correspond to the 1% significance level of variability evaluated accounting for the S/N of the data. On the sides of the Mg II line, significant narrow variabilities are due to the imperfect correction of the telluric lines.
In the text

Fig. 7. Comparisons between the centroid of He I λ 6678 Å and the centroids of other lines detected in the spectrum of Schulte 12 (He I lines on top left, H I lines on top right, and metallic absorption and emission lines on bottom left and right, respectively). Filled symbols correspond to Hermes data, open symbols and crosses to Carmenes data; the diagonal yields the one-to-one correlation. For comparison, the case of a DIB is added to the last panel.

In the text

	Fig. 8. Evolution of the photometry of Schulte 12 and of the centroids of some prominent lines in its spectrum during the Carmenes campaign in Summer 2018.
In the text

	Fig. 9. Periodograms calculated on the line profiles (see text) for H I λλ 8598, 8663, 10 049, 15 556 Å, He I λλ 5876, 6678, 7065 Å, C II λλ 6578, 6583 Å, O I λλ 7772, 7774 Å, Mg II λ 10 952 Å, and the DIB at 6613 Å for comparison. The vertical dotted line shows the frequency corresponding to the 108 d period.
In the text

Fig. 10. Tentative orbital solution of Schulte 12 for a fixed orbital period of 107.99 d (middle column of Table 4). The average velocities of Mg II and Fe II emissions in the Carmenes dataset are shown by blue dots and the Hermes velocities of Fe II λ 7513 Å (not used to obtain the orbital solution) by open red squares. We note that the phase used here corresponds to the binary ephemeris (Table 4), not to the X-ray ephemeris used for Figs. 2 and 3.

In the text

Fig. 11. Revised position of Schulte 12 in the Hertzsprung-Russell diagram adopting the distance corresponding to the Gaia-DR2 parallax. The evolutionary tracks are taken from Ekström et al. (2012) for solar metallicity and no rotation. The continuous red line corresponds to the empirical Humphreys–Davidson limit (Humphreys & Davidson 1994). Various instability boundaries are shown according to Saio (2011). The dotted red curve indicates the instability boundary of low-order radial and non-radial p-mode oscillations. The short-dashed magenta contour indicates the region where non-radial g-modes with ℓ ≤ 2 are excited by the iron opacity bump. Finally, the cyan long-dashed contour indicates the region where oscillatory convection modes with ℓ ≤ 2 are excited (Saio 2011).

In the text

	Fig. A.1. X-ray spectrum of Schulte 12 (left: Swift data of observation 00032767003, right: XMM-Newton data taken on rev. 3176 with pn in green, MOS1 in black, and MOS2 in red) along with their best-fit models (see Table 2) and the residuals.
In the text

	Fig. A.2. Results from period searches on the visible light curves. From left to right: results for OMC data, ASAS-SN data (all runs and last run only), and our own photometry. In each panel, from top to bottom: periodogram resulting from the Fourier method (with its spectral window), the AOV method, and the conditional entropy method (see text for details).
In the text

	Fig. A.3. Variation of the line strengths and the photometry in 2018. Formal errors on EWs are similar to the symbol size. The dashed cyan and dotted black vertical lines indicate times of maximum and minimum brightness, respectively.
In the text