Asteroseismology of Massive Stars with the TESS Mission: The Runaway β Cep Pulsator PHL 346 = HN Aqr

PHL 346 is a star of V = 11.44 located at a Galactic latitude of b ≈ 58°. Kilkenny et al. (1977) classified it as spectral type B1, and Keenan et al. (1986) reported a surface gravity that is consistent with an evolved main-sequence evolutionary status and Pop. I metal abundances. Given the radial velocity that they measured (+66 ± 10 km s⁻¹), Keenan et al. (1986) had to conclude that PHL 346 would not have had enough time to attain such a high Galactic latitude within its lifetime, and thus must have been formed far from the Galactic plane. This puzzling result was amended by Ramspeck et al. (2001) who, based on a new spectroscopic analysis and the first proper motion measurement of PHL 346, reconciled the stellar flight time with its lifetime, meaning that the star could have been formed in and ejected from the Galactic plane.

PHL 346 was the subject of several studies that resulted in determinations of its effective temperature and surface gravity, summarized in Table 2 and available online only. β Cep-type pulsations of PHL 346 were discovered by Waelkens & Rufener (1988) and confirmed by Kilkenny & van Wyk (1990). The star was also observed during the ASAS survey (Pigulski & Pojmański 2008) and by Handler & Shobbrook (2008). All these authors detected the same single oscillation frequency near 6.566 day⁻¹. Heynderickx et al. (1994) and Cugier et al. (1994) suggested that this pulsation is due to a nonradial l = 1 mode. On the other hand, Handler & Shobbrook (2008) derived l = 2 or 4 for this oscillation, but noted a possible problem with the U filter data that their identification critically hinged upon.

Given that HN Aqr is the first β Cep pulsator observed with TESS, that its light variations are multiperiodic, and that the time base of the observations may not allow all pulsational signals to be resolved in frequency, we analyzed the data with various methods. Method 1 comprised classical single-frequency power spectrum analysis and simultaneous multifrequency sine-wave fitting with input parameter optimization. The sine-wave fits are subtracted from the data and the residuals are examined for the presence of further periodicities ("prewhitening"). The noise level was assessed in apparently signal-free frequency regions and the signal-to-noise ratio (S/N) >4 criterion by Breger et al. (1993) was adopted to evaluate the significance of signal detection.

Method 2 used essentially the same approach, but the background noise in the Fourier spectrum was modeled in log–log space following Pablo et al. (2017) using the amplitude spectrum from 0.036 to 30 day⁻¹ prewhitened with the strongest oscillation. Again, the S/N > 4 criterion was used. Method 3 applied the MIARMA gap-filling method (Pascual-Granado et al. 2015) to improve the spectral window function of the data. This modified light curve was then frequency analyzed with the SigSpec algorithm (Reegen 2007) that also involves prewhitening.

Method 4 (Kallinger & Weiss 2017) employed a Bayesian algorithm allowing a probabilistic assessment of the significance of signal detection. It was run twice, once allowing only a single frequency to represent a peak in the amplitude spectrum (equivalent to multifrequency sine-wave fitting), and the other time allowing close frequency doublets (spacing ≤3/ΔT = 0.11 day⁻¹). Finally, Method 5 used a Morley wavelet transform (Torrence & Compo 1998). Output from this approach reflected the amplitude modulation of the short-period pulsations (Figure 1). In the low-frequency region (f < 2 day⁻¹), considerable changes in the amplitudes of individual signals were visible, in most cases with some degree of regularity indicating multifrequency beating.

To reconcile the outcome of the different methods we have compared the results of the first four techniques that resulted in lists of frequencies and amplitudes. They yielded fairly consistent results for strong, well-resolved signals, but diverged in frequency regions where densely spaced signals were present. The number of detected frequencies strongly depended on the adopted signal detection threshold.

To provide a set of pulsation frequencies that can be reasonably safely used for asteroseismic investigations, we only accepted signals detected by at least three of the methods independently. The amplitudes of these signals had to exceed the noise level a(ν) (Equation (1)) by factors of 5 (independent signals) and 3.5 (combination frequencies), respectively. The more conservative S/N threshold was chosen to avoid picking up spurious frequencies in a data set of the present size and sampling (Baran et al. 2015); the noise level was computed according to Method 2:

$$\begin{eqnarray}&&a(\nu )=\displaystyle \frac{{a}_{0}}{1+{\left(\tfrac{\nu }{{\nu }_{0}}\right)}^{\lambda }}+{P}_{0},\end{eqnarray} \tag{ 1 }$$

where P₀ = 1.4388 × 10⁻⁵ mag is a constant Gaussian noise term, λ = 2.0, a₀ = 7.76 × 10⁻⁵ mag and ${\nu }_{0}={(2\pi \tau )}^{-1}$, τ = 0.076 day are the amplitude and characteristic frequency describing the red noise, respectively.

The resulting list of frequencies is given in Table 1. The frequency values and error estimates were adopted from the Bayesian method as it does not rely on prewhitening. The error bars correspond well to the least-squares errors (Montgomery & O'Donoghue 1999) for well-separated signals, but take systematic uncertainties of closely spaced frequencies into account. As it is not yet well known how much data processing affects the low-frequency domain, periods longer than 4 days should be treated with caution. Some steps of prewhitening of the TESS data are illustrated in Figure 2. Our multifrequency fit leaves a residual scatter of 3.4 mmag per point (617 ppm hr⁻¹) containing residual systematic (presumably mostly stellar) variability contributing some 15% to the total scatter.

Zoom In Zoom Out Reset image size

The low-frequency domain of detected signals extends up to 1.49 day⁻¹. As the large number of frequencies cannot be explained by effects of binarity or rotation most, if not all, of these must be due to pulsation. There is no evidence of rotational or binary-induced modulation, or of a regular period spacing. A re-occurring frequency spacing is present within the signals in the p-mode domain (Δf ≈ 0.187 day⁻¹), which may be due to rotational splitting. No multiplets were identified with confidence; the best candidate consists of frequencies f₂₂ − f₂₅, which could be part of an l = 2 quintuplet with the m = 0 mode missing.

Thirty years have passed since pulsation in PHL 346 was discovered. Thus the long-term stability of the period and amplitude of the dominant pulsation mode can be studied. We gathered all available time-series photometry of the star (Waelkens & Rufener 1988; Kilkenny & van Wyk 1990; Handler & Shobbrook 2008; Pigulski & Pojmański 2008). Additionally, we include some previously unpublished observations from 1989 and 1990 by one of us (D. K.), International Ultraviolet Explorer (IUE) spectrophotometry (Dufton et al. 1998), as well as Northern Sky Variability Survey (NSVS; Woźniak et al. 2004), and All Sky Automated Survey for SuperNovae (ASAS-SN; Shappee et al. 2014) data. Each of the individual data sets used spans at least 13 days.

The resulting O–C diagram and amplitudes over time are shown in Figure 3. Taking into account passband differences between the various individual data sets, the amplitude of the dominant mode remained relatively stable over 30 yr at a level of about 22 mmag and dropped only in the late 1980s and early 1990s, down to 11.9 ± 0.7 mmag in 1990. (As expected, the amplitude in the ultraviolet (UV) was higher than in the visual.) During this drop, the O–C diagram shows a significant period change: the 1990 value of O–C is ≈25 minutes off the ephemeris. We checked whether the shape of the O–C diagram can be explained in terms of the light-time effect in a binary system, with no satisfactory solution. In any case, such an orbit needs to be highly eccentric (e ≳ 0.7) with periastron passage in the mid 1990s, when unfortunately a gap of almost a decade in the photometric data occurs. This would be also the time when large radial-velocity (RV) changes would occur. We therefore gathered all available RV data for HN Aqr (Keenan et al. 1986; Kilkenny & Muller 1989; Hambly et al. 1996; Ramspeck et al. 2001; Lynn et al. 2002), and plot them in the lower panel of Figure 3. Although there is no large change of RV at the predicted time of the periastron passage, the first two RVs (Kilkenny & Muller 1989), are significantly different.

Zoom In Zoom Out Reset image size

Assuming a period of 30 yr for a hypothetical binary orbit of HN Aqr and M = 10 M_⊙ for the primary, a companion of M₂ sin i = 1.15 M_⊙ would be able to explain the 0.013 day semi-amplitude of the (O–C) variations. However, such a companion would cause an orbital RV semi-amplitude of only 2.2 km s⁻¹, more than an order of magnitude less than the observed RV change. To summarize, in view of the large uncertainties of RVs, the presence of pulsations that increase the scatter of RVs, the large gap in photometric data that causes cycle count ambiguities in the 1990s, the inconsistency of the hypothetical light time effect, and the radial velocity change, the results as to the binarity of the star are inconclusive. What one can conclude, however, is that the observed pulsational (O–C) variations and RV changes cannot be caused by binarity alone.

High-resolution (R ≈ 35,000) spectra of HN Aqr were taken with the Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) attached to the Very Large Telescope (VLT)-UT2 at the European Southern Observatory on the night of 2009 August 27. Integration times of 1150 s at both the blue and the red arm led to S/N of about 150–190 in the range of 304–669 nm, and around 70 up to 1043 nm. The parameters derived from these spectra using two analysis strategies (Irrgang et al. 2014; Hubeny & Lanz 2017) are listed in online Table 2.

Summarizing all these results and using bolometric corrections from Flower (1996) suggests that a range of basic parameters T_eff = 22,300 ± 900 K, log g = 3.75 ± 0.15, M_bol = −5.2 ± 0.3 and [M/H] = 0.2 ± 0.1 should contain the parameter space in which a seismic model for HN Aqr ought to be located. Furthermore, a projected rotational velocity v sin i = 30 ± 5 km s⁻¹ can be adopted. The Gaia second data release (DR2) parallax π = 0.10 ± 0.08 mas for HN Aqr (Gaia Collaboration et al. 2016, 2018; Luri et al. 2018) cannot be used to better constrain the stellar luminosity, but is consistent with the results above.

The increased metallicity agrees very well with a kinematic investigation based on the radial velocity, Gaia DR2 proper motions, and a derived spectrophotometric distance of 6.7 ± 0.8 kpc (99% confidence interval; assuming a mass of 9.5 ± 0.4 M_⊙). It suggests that the star stems from the inner part of the Galaxy (Galactocentric radius at disk intersection 2.3 ± 0.4 kpc with 68% confidence interval). After taking the Galactic abundance gradients into consideration, the abundance pattern appears to be perfectly normal except for an underabundance of ∼0.3 dex in carbon. From the calculation of stellar trajectories (for details, see Irrgang et al. 2013), we infer a flight time to the Galactic disk of 23 ± 1 Myr (68% confidence interval).

We follow Daszyńska-Daszkiewicz et al. (2017) to model the pulsation spectrum of HN Aqr preliminarily. The basic stellar parameters determined earlier are represented best by models of 11 M_⊙ that we take for the purpose of example. We attempt to reproduce the frequency domains of the observed p- and g-modes (Section 3.1) and search for required modifications of the stellar models. Given the low projected rotational velocity of the star, and that the suspected first-order p-mode splitting from Section 3.1 would suggest v_rot ≈ 64 km s⁻¹, the observed frequency ranges should not differ much from those in the stellar frame of rest. In addition to the frequencies, we consider the age constraint of 23 ± 1 Myr from the kinematic analysis. The most important influence on mode stability is provided by the overall metallicity, as shown in Figure 4.

Zoom In Zoom Out Reset image size

However, increasing Z alone, consistent with the spectroscopically determined abundances (Table 2), which gives a better match to the p-mode domain, is insufficient to reproduce the observed g-mode region. Also, increasing Z requires more massive models that evolve too rapidly: the model with Z = 0.025 has an age of 13.5 Myr and a main sequence lifetime of 18.0 Myr. Possible ways out of this problem would be the inclusion of convective core overshooting prolonging the main sequence lifetime (a 11 M_⊙, Z = 0.025, α_ov = 0.2, $\mathrm{log}{T}_{\mathrm{eff}}=4.35$ model has an age of 15.5 Myr), increasing the He abundance, leading to models of lower mass, and modifications of the input opacities, in particular near $\mathrm{log}T=5.46$, corresponding to an enhancement of the nickel opacity (see, e.g., Daszyńska-Daszkiewicz et al. 2017).