© ESO, 2015

1. Introduction

The source 2M 1938+4603 (J19385+4603, KIC 9472174) was recognized as a binary system consisting of a subdwarf B (sdB) star and an M dwarf (Østensen et al. 2010). The sdB star shows pulsations, and it is potentially a prominent target for asteroseismology. The light curve of the system exhibits two shallow minima as a result of mutual eclipses of the components of the binary system and a very strong reflection effect. With two kinds of variations, geometric and pulsational, it might be possible to use both in combination to better constrain the physical properties’ primary component, which makes 2M 1938+4603 a promising system for asteroseismology purposes.

Up to now, a dozen systems similar to 2M 1938+4603 have been found (Zorotovic & Schreiber 2013). The orbital period of 3 h falls in the range of the other system’s periods. However, the difference from the previously known systems is seen in an amplitude of the reflection effect and grazing eclipses that dominate the binary variation. These systems are usually the object of light curve analysis as well as O−C analysis when searching for the presence of low-mass tertiary components.

The results obtained from analyzing the Q0 data along with earlier ground observations were published by Østensen et al. (2010). The authors made an attempt to model the orbital trend with the ultimate goal of analyzing pulsations tied to the sdB star. It turned out, however, that the reflection effect is too complex to be modeled with the current approaches, so they decided to mathematically remove the trend, instead. They concluded that more sophisticated atmosphere models are needed to better describe the heating of the secondary, which is responsible for such a strong reflection effect. The authors constrained the masses of the components from the spectroscopic observations of the radial velocity amplitude and the surface gravity of the sdB star. A reliable value for the inclination angle was also found from modeling of the eclipses in the light curve. After removing the binary trend, Østensen et al. (2010) analyzed the residual data and obtained a list of pulsation frequencies for future asteroseismic analysis. They listed 55 frequencies that were significant in the nine-day Q0 data; however, some of them are likely combination frequencies and orbital aliases. The authors also used midtimes of 13 primary eclipses obtained from the ground over a two-year period to derive a linear ephemeris.

Soon after, Barlow et al. (2012) used a longer data coverage (Q0,5,6) to constrain masses of the components by means of the Rømer delay. The authors derived a value for that delay of around 2 s. The masses they derived were different from those obtained by Østensen et al. (2010). However, the values from these two analyses could be aligned by adjusting other parameters, since the modeling of the binary trend is usually affected by degeneracies among the parameters involved in the solutions. The authors used midtimes of eclipses to construct the O−C diagram. They removed the pulsation signal, using the list of frequencies given by Østensen et al. (2010). These low amplitude variations affect the midtimes in a significant way, giving rise to many spurious periods in the O−C diagram, with ~3.7 days being the most significant. According to their analysis, there is an indication of evolutionary orbital period decay. The authors notice that this decay is several orders of magnitude faster than the values observed in other similar systems.

In this analysis, we used the near-continuous photometry of 2M 1938+4603 collected during the first stage of the Kepler mission. While we initially intended to investigate the pulsation properties of the sdB star to describe it fully in similar terms to the many other interesting pulsators observed by the Kepler spacecraft (e.g., Baran et al. 2012; Telting et al. 2012; Reed et al. 2014; Østensen et al. 2014) and hopefully to be able to interpret the pulsation spectrum in terms of rotational multiplets, this was not feasible because of the complexity of the period and amplitude variability in many of the modes. For the purpose of deriving the eclipse timings, we were satisfied with removing only the most significant modes in the pulsation spectrum so that reliable eclipse timings could be measured and analyzed.

2. Data

The source 2M 1938+4603 was observed during Q0 and Q5-17 in the short cadence (SC) mode, as well as Q0-17 in the long cadence (LC) mode. We downloaded all data from the Barbara A. Mikulski Archive for Space Telescopes (MAST)¹. We worked both with pixel data that gives access to all individual pixels within the target mask and with simple aperture photometry (SAP) data. The latter represents the sum of the fluxes from all pixels within the optimal aperture for the current target. We have not noticed any significant contamination and have accepted that the SAP fluxes are representative of the best approach, since no other pixel combination gave a better signal-to-noise ratio.

It is commonly known that the data obtained with the Kepler spacecraft are affected by many systematics that have to be removed prior to the final analysis. For this purpose we used the multitask PyKE package, particularly the kepcotrend task with a set of cotrending basis vectors provided by the mission specialists. The additional long-term detrending we applied prior to the Fourier analysis is described in Sect. 3.1. Next, we clipped data at 4σ and converted them to parts per notation. Finally, all individual months were stitched together for further analysis.

Fig. 1 Light curve of KIC 9472174 phased over the orbital period with pulsations removed.

3. Flux variations

The light curve of 2M 1938+4603 has two types of variations: geometrical light changes caused by mutual orbiting of the components and a physical variation as a result of stellar oscillations present in the sdB star. In principle, these two types of variations can provide an independent estimation of the physical properties of the hotter component. However, it is very difficult to analyze them simultaneously, so they have to be separated before detailed analysis of each type. Since the binary dominates the light-curve variations, we first attempted to remove it from the data in order to leave a clean pulsation signal. We present the plot showing the binary trend with pulsations removed in Fig. 1.

3.1. Binary solution

The code most commonly used for light curve modeling is the Wilson-Devinney (W-D) model (Wilson & Devinney 1971). However, this code does not work well in the case of a cool star strongly heated by a very hot companion, even if the enhanced reflection model is being used. At present, a grid of models of the heated atmosphere is not available to us. Therefore, we only made some alterations to the original W-D model by using the limb darkening coefficients for the Kepler passband (from Claret et al. 2013) according to the local temperature of the heated hemisphere of the secondary star, which was significantly increased by the hot companion radiation. We also treated the albedo of the secondary star as a free parameter and allowed it to be greater than one. These two alterations yielded a pretty good description of the overall shape of the light curve; however, some resulting parameters turned out to be unphysical: we obtained 16 000 K for the effective temperature of the sdB star and about 2 for the albedo of the secondary. A more detailed explanation of the binary modeling of 2M 1938+4603 is provided in Zola & Baran (2013).

Fig. 2 Amplitude spectrum calculated from rectified Q5-Q17 data limited to the frequency range with signal detected. The horizontal dashed line denotes our arbitrarily adopted threshold for frequency removal.

Since the models we came up with are unphysical, we could not use them to derive the physical properties of the system. Instead, we used these best models, computed for each month separately, only as representations of the binary variation to remove long-term trends that were not removed by the kepcotrend procedure (Kinemuchi et al. 2012). To do so, we calculated differences between actual data and the models, then we fitted a cubic spline curve to these residuals and removed that fit from the original (not residual) data. This detrending removed a systematic scatter from phased data, which were used to remove the binary trend. We did it in the following way. We phased data with a preliminary ephemeris and fitted a cubic spline curve to these data. Then, this fit was subtracted from the unphased data. The remaining data were free of systematics and the binary trend and ready for the Fourier analysis.

3.2. Oscillations

The data prepared in the way described in the previous section were subject to Fourier analysis. We present the amplitude spectrum in Fig. 2. The formal frequency resolution, defined by 1.5 /T, equals 0.026 μHz, while the mean noise level N = 0.5 ppm. The amplitude threshold for removing frequencies that significantly contaminate the signal with respect to the following O−C analysis was adopted at 0.011 ppt.

We searched the amplitude spectrum for rotationally split multiplets and sequences evenly spaced in period, which could help for understanding the pulsation spectrum observed in the primary component. We found no multiplets, while the Kolmogorov-Smirnov test, commonly used in the analysis of other sdBV stars observed with the Kepler spacecraft (e.g., Reed et al. 2011), revealed no indication of modes evenly spaced in period. Unluckily, the spectrum is very complex, because it contains fake peaks arising from flux modulation, and we were unable to fully analyze the flux variation caused by pulsations.

4. Stability of the orbital period

We used the eclipses observed in the light curve to study the stability of the orbital period of the binary system. A mid-time of each minimum can be compared with the one calculated assuming a constant period and plotted along time. The resulting diagram, often called observed-minus-calculated (O−C), is a very handy tool for verifying period stability (further information on the O−C analysis can be found in Sterken 2005).

As already shown by Barlow et al. (2012), the pulsations present in the primary component distort the shapes of minima significantly, leading to improper estimation of midtimes. By adding pulsations to simulated light curves, Barlow et al. (2012) demonstrated that many peaks between 3.6 days and 0.3 days are caused by the pulsations and are not true orbital effects.

Fig. 3 Top panel: the O−C diagram for primary and secondary minima based on mid-points derived with the KW method (red − primary minima, black − secondary minima). Bottom panel: a comparison of primary minima derived with a variety of methods (red = KW, blue = parabola, green = model). See Sect. 4 for more details.

We used prewhitening to remove 58 pulsation periods with amplitudes higher than 0.011 ppt from the data. We considered the remaining flux variations having negligible effect on the midtimes derivations. In fact, the spurious 3.7-day period did not show up in the O−C diagram. We then used the method described in Sect.  2 of Kwee & van Woerden (1956) (hereafter: KW) to derive the eclipse mid-point times. We assigned each mid-point to a proper epoch and fitted a linear regression to estimate the average orbital period. The linear ephemeris is ${\mathit{T}}_{\mathrm{prim}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{2}\hspace{0.17em}\mathrm{455}\hspace{0.17em}\mathrm{276.60843}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathrm{+}\mathit{E}\mathrm{·}\mathrm{0.125765282}\mathrm{\left(}\mathrm{5}\mathrm{\right)}\hspace{0.17em}\mathrm{d}$(1)where E represents consecutive epochs, while prim stands for primary.

Fig. 4 Top panel: the O−C diagram for primary minima based on mid-points derived with the KW method overplotted with the three fits we applied. The blue line represents a parabola+sine fit, the red line represents a sine+sine+line fit, and the green one shows a sine+sine fit. Bottom panel: residuals after the sine+sine+line fit removal.

In the top panel of Fig. 3 we present the O−C diagram of the primary and secondary minima obtained by means of the KW method. Although the secondary minima are noisier than the primary ones, it is clear that the long-term variations follow the same trend, which excludes an apsidal motion as the cause of the O−C variation. To reject systematics as a potential source of that variation, we calculated the O−C with alternative methods, by applying the KW method to the mean of ten consecutive midtimes, fitting parabolae to the minima and fitting the phase shift parameter (defined in Wilson & Devinney 1971) in monthly sets of data used for the binary modeling by Zola & Baran (2013). The O−C diagrams obtained with a variety of methods are plotted in the bottom panel of Fig. 3. As can be seen, all diagrams provided consistent results, though with different levels of precision, therefore we assume that the O−C variations are intrinsic to the binary system.

Sterken (2005) lists potential shapes of variations in an O−C diagram and explained their astrophysical contexts. A shape can be linear, which is indicative of a constant period that has still not been correctly adopted. It can also be parabolic (or irregular) as a consequence of evolutionary changes in the system or sine-like, if caused by a companion(s) to a binary system. A combination of many shapes may cause complexity in the interpretation of possible sources of the O−C variation. In the case of 2M 1938+4603, we at first interpreted the O−C variation as a mix of a parabola and a sinusoid. The parabola was open upward, an indication that the orbital period is increasing. We had no explanation of the source of this increase and it could also be that the parabola was, in fact, a part of another long-term sinusoidal variation. A similar behavior of an increasing period was found in HS 0705+6700 (Qian et al. 2013), and since neither gravitational waves radiation nor magnetic braking can lead to this, the authors speculated about another companion to that system.

Then we fitted two sinusoidal components (with and without a linear trend to account for an imperfect period estimation). Since two sinusoids with a linear term fit reproduces the shape of the O−C well, and we have no explanation for an open-up parabola, we followed Qian et al. (2013) and decided to accept that fit as the best solution. It also turned out that this fit (among all three we used) predicts Q0 data, which were not included in our analysis so far, very well. In the final solution we included all the data available to us, including Q0. To avoid any unnecessary influence of one trend on another, we fitted both sinusoids simultaneously. We only used the primary minima, which are characterized by lower noise. We stress that the longer-period sinusoid should be considered with caution. Its amplitude and period increased significantly, which is a consequence of the fact that not even one cycle of this variation is covered. More data will likely change this sine term even further, therefore, it is too early to draw any astrophysical conclusions based on the second sinusoid. We present a comparison of all three fits in the top panel of Fig. 4 while residuals after the best fit removal are added to the bottom panel.

As an indication of a third body orbiting two inner components, the sinusoid has an amplitude of 1.27(2) s and a period of 416(2) days. The longer-term sinusoid has an amplitude of 3.2(2) s and a period of 1797(51) days. The uncertainties represent standard 1σ fitting errors.

To estimate the mass of the tertiary companion, we need to know the masses of the two stars in the binary system. To derive these two values we made an effort to estimate the Rømer delay in the same manner as in Barlow et al. (2012).

5. The offset of the secondary minimum

Kaplan (2010) shows that the finite speed of light can be used to determine the physical size of the orbit and to constrain masses of the components in the binary system NLTT 11748. In the systems with two components of unequal masses, the light travel delay (aka Rømer delay) causes the secondary minimum to occur at not exactly half the orbital period. According to Eq. (1) in Kaplan (2010), the size of this delay is a direct measure of the semi major axis and the total mass of the system. Barlow et al. (2012) made the first attempt to use this feature to estimate the masses of the two components in 2M 1938+4603. They derived the delay to be around 2 s and subsequently the mass ratio of the system as 0.2691.

We determined an offset of the secondary minimum for each timestamp. We did not account for the orbital period variation since the period change during one orbital revolution is negligible and can affect the Rømer delay by less than 3 × 10^-4 s. Based on the longer time span, the average Rømer delay equals 1.70(7) s, where the uncertainty represents 1σ fitting error. It should be mentioned here that when using the same time-span as available to Barlow et al. (2012), we obtained the same result as they did.

Having determined the average light travel delay, the binary period and the amplitude of the radial velocity of the sdB star, the mass ratio of the system can be evaluated. It has been already done by Barlow et al. (2012), although the mass of the primary component of 0.3 M_⊙ is rather low for sdB stars. We even derived a lower Rømer delay value that implies that the sdB mass would be lower than 0.3 M_⊙. Such a value does not agree with the gravity estimated from spectroscopy (Østensen et al. 2010). This leads to the paradox that more data gives a worse estimation of the sdB mass, therefore we started to question the correctness of the Rømer delay estimation in this case. Barlow et al. (2012) have already mentioned that both the Rømer effect and the eccentricity (e) shift the secondary minimum off the phase 0.5 and a very small e would reproduce the observed shift in the secondary minimum. As such, the two effects become indistinguishable. However, while the Rømer effect keeps the shape of the minimum symmetric, an elliptical orbit will cause its distortion. Unfortunately, owing to a very small eccentricity, that distortion will be too small to be seen in the light curve. In addition, while the Rømer delay always shifts the mid-time of the minimum to a larger phase, the eccentricity can shift it either way. This means that the overall shift of the minimum can be a superposition of the two effects, yet their disentangling is hardly possible. To confirm our suspicion, we played with an eccentricity to derive the same offset of the secondary minima as associated with the Rømer delay. We derived e = 0.0004, which is comparable to the predictions made by Barlow et al. (2012).

To reach 0.48 M_⊙ for the sdB star and q = 0.25 (after Østensen et al. 2010), the Rømer delay needs to be 2.27 s. We expect the Rømer delay to be present in this system because the masses of the two main components are not equal. Since we derived a value of 1.7 s, the only way to have a mass of the sdB star in agreement with the spectroscopic gravity determination is to postulate that a non-zero eccentricity must also be present, shifting the secondary minimum in the opposite direction. Using Eq. (2) from Barlow et al. (2012) we estimated ecosω = 0.00008, where ω is poorly determined for small eccentricities, though the value of cosω must be negative. We stress that in case of such small eccentricity, it is impossible to measure it from the currently available data (both photometric and spectroscopic).

Our analysis clearly shows that a reliable determination of the Rømer delay in the available Kepler data of 2M 1938+4603 is currently impossible. The derived value of this delay, without accounting for a non-zero eccentricity, leads to unreliabling estimating the properties of this binary system. Therefore, the mass of the tertiary component can only be derived by adopting the necessary parameters from Østensen et al. (2010). We adopted the masses of the components M₁ = 0.48 M_⊙ for the sdB star and M₂ = 0.12 M_⊙ for the M dwarf. We employed the third Kepler law to estimate the size of the circular planetary orbit from the barycenter of the binary stars, deriving 0.92(2) AU. Next, we used the amplitude of the O−C variation of 1.30 s and derived the minimum planetary mass M₃·sini = 1.8(1) Jovian masses. Assuming that the planetary orbit is coplanar with the binary orbit, the planetary mass will reach 1.9(1) Jovian masses. The mass of a brown dwarf can only be reached if the inclination of the planet’s orbit deviates substantially from the plane of the inner binary and is as low as 7.1 deg. We provide the parameters of the entire triple system under the coplanar assumption in Table 1.

We conclude that 2M 1938+4603 is an sdB+M system with a circumbinary companion. A brown dwarf tertiary was found in HS 0705+6700 (Qian et al. 2009, 2012), while planets were proposed to orbit, for example, HW Vir (Kilkenny et al. 2003; Lee et al. 2009; Beuermann et al. 2012), NSVS 14256825 (Almeida et al. 2013; Hinse et al. 2014), and NY Vir (Qian et al. 2012).

6. Summary

We have presented our analysis of Q5-17 Kepler photometric data of the triple system 2M 1938+4603, consisting of an sdB star, an M dwarf, and a planet orbiting the binary. The light curve shows both orbital effect from the inner binary and pulsations from the sdB primary. We could not analyze both effects simultaneously, so we attempted to treat them separately. First we modeled the binary variations, as described in Zola & Baran (2013). As can be seen, the shape of the light curve is reproduced fairly well, but the final model is unphysical, so the resulting system parameters cannot be relied on. We made use of this model as a first-order representation of the orbital light-curve variations to remove them before analyzing the pulsations.

We calculated the amplitude spectrum from data spanning more than 37 months, hoping to be able to make some asteroseismic mode identification. We looked for multiplets as a consequence of stellar rotation and gravity modes evenly spaced in periods. Unfortunately, we found none. The lack of multiplets has been already found in another sdB star observed with Kepler (Baran et al. 2015). The Kolmogorov-Smirnov test gave no clear evidence of gravity mode sequences that are evenly spaced in period. In addition, there were problems with doing reliable prewhitening, most likely because of amplitude/phase variability. A complete list of peaks was therefore unfeasible, so we decided to focus solely on the eclipse timing variations.

We used the midtimes of the primary minima to verify the stability of the orbital period of the system. We removed the highest amplitude pulsations and split the data into individual minima. The O−C diagram of the orbital period clearly shows that the period is not constant. We fitted both a parabolic and a sinusoidal trend to the data. The former is often associated with an evolutionary decrease in the period, while the latter may indicate an additional companion to the system. We derived an open-up parabola and had no explanation for this variation, except that it might be part of a very long sinusoidal variation. The sinusoid was interpreted as being due to a planet orbiting the system. Our discovery yields that 2M 1938+4603 is another sdB+M system with a substellar companion and one of the lowest masses among all tertiaries in other well-studied sdB+M binaries.

Finally, we used the phase shift between primary and secondary eclipses to estimate the Rømer delay. This effect allows measuring the masses of the components. While this analysis is valid for a circular orbit, the resulting low mass is unconvincing as a very small eccentricity could easily reduce the Rømer delay by the required amount. Therefore, we can only conclude that the inner binary can either be a low-mass subdwarf with a substellar companion in a circular orbit or a normal sdB+dM binary in a slightly elliptical orbit. If the latter is the case, continued observations of the Rømer delay over the coming decades should reveal a change in the delay as the orientation of the ellipse precesses. But this may be very hard to achieve without reliable disentanglement of the pulsations from the eclipsing light curve. We hope that the discovery of a Jupiter-mass planet orbiting the system encourages people to continue observing this system from the ground.

Acknowledgments

This work was supported by National Science Center under projects No. UMO-2011/03/D/ST9/01914 and 2012/07/B/ST9/04432. A.S.B. thanks Wojtek Dziembowski for reading a draft and for his useful comments. We thank the entire Kepler team for the many years of work that is proving so successful.

References

All Tables

All Figures

	Fig. 1 Light curve of KIC 9472174 phased over the orbital period with pulsations removed.
In the text

	Fig. 2 Amplitude spectrum calculated from rectified Q5-Q17 data limited to the frequency range with signal detected. The horizontal dashed line denotes our arbitrarily adopted threshold for frequency removal.
In the text

	Fig. 3 Top panel: the O−C diagram for primary and secondary minima based on mid-points derived with the KW method (red − primary minima, black − secondary minima). Bottom panel: a comparison of primary minima derived with a variety of methods (red = KW, blue = parabola, green = model). See Sect. 4 for more details.
In the text

	Fig. 4 Top panel: the O−C diagram for primary minima based on mid-points derived with the KW method overplotted with the three fits we applied. The blue line represents a parabola+sine fit, the red line represents a sine+sine+line fit, and the green one shows a sine+sine fit. Bottom panel: residuals after the sine+sine+line fit removal.
In the text