The EXOTIME project: Signals in the $ O-C $ diagrams of the rapidly pulsating subdwarfs DW Lyn, V1636 Ori, QQ Vir, and V541 Hya

We aim to investigate variations in the arrival time of coherent stellar pulsations due to the light-travel time effect to test for the presence of sub-stellar companions. Those companions are the key to one possible formation scenario of apparently single sub-dwarf B stars. We made use of an extensive set of ground-based observations of the four large amplitude p-mode pulsators DW Lyn, V1636 Ori, QQ Vir, and V541 Hya. Observations of the TESS space telescope are available on two of the targets. The timing method compares the phase of sinusoidal fits to the full multi-epoch light curves with phases from the fit of a number of subsets of the original time series. Observations of the TESS mission do not sample the pulsations well enough to be useful due to the (currently) fixed two-minute cadence. From the ground-based observations, we infer evolutionary parameters from the arrival times. The residual signals show many statistically significant periodic signals, but no clear evidence for changes in arrival time induced by sub-stellar companions. The signals can be explained partly by mode beating effects. We derive upper limits on companion masses set by the observational campaign.


Introduction
Subdwarf B stars (sdBs) are sub-luminous stars with a mass of about 0.5 M located at the blue end of the horizontal branch, which is the so-called extreme horizontal branch (EHB, Heber Based on observations obtained at the 0.9 m SARA-KP telescope, which is operated by the Southeastern Association for Research in Astronomy (saraobservatory.org).
Photometric data of Fig. 1,results in Figs. 8,10,12,and 15, and figures in the appendix are available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via . e-mail: mackebrandt@mps.mpg.de 1986). They maintain a helium burning core, but their thin hydrogen envelope (M env < 0.01 M ) cannot sustain hydrogen shell burning, which identifies sdBs with stripped cores of red giants (Heber 2016). Binary evolution with a common envelope (CE) is the favoured formation scenario for most sdBs. The sdB progenitor fills its Roche lobe near the tip of the RGB. A CE is formed when the mass transfer rate is sufficiently high and the companion star cannot accrete all the matter. For close binary systems with small initial mass ratios q < 1.2-1.5, two phases of mass transfer occur. The first Roche-lobe overflow is stable, whereas the second one is unstable, leading to the ejection of the CE. The resulting binary consists of an sdB star and a white dwarf in a short-period orbit. For initial mass ratios q > 1.2-1.5, the first mass-transfer phase is unstable and the CE is ejected, producing an sdB star with a non-degenerate (e.g. main sequence star) companion. A more detailed review of formation and evolution of compact binary systems can be found in Podsiadlowski (2008) and Postnov & Yungelson (2014). These formation scenarios cannot explain the observed additional occurrence of apparently single sdBs (Maxted et al. 2001). Among the proposed formation scenarios is the proposal by Webbink (1984) that they could be formed by a merger of two helium white dwarfs. But such mergers are problematic, as they are expected to retain very little hydrogen (Han et al. 2002) and be left with higher rotation rates than what has been observed . Moreover, the overall observed mass distribution of single sdB stars is not consistent with that expected from the proposed formation scenario. Sub-stellar companions could resolve this disagreement between theory and observations. Planetary-mass companions like the candidates V391 Peg b (Silvotti et al. 2007), KIC 05807616 b,c (Charpinet et al. 2011), KIC 10001893 b,c,d (Silvotti et al. 2014), or brown dwarf companions like V2008-1753 B (Schaffenroth et al. 2015) or CS 1246 (Barlow et al. 2011b) indicate the existence of a previously undiscovered population of companions to apparently single sdBs. Due to the high surface gravity and effective temperature (leading to few, strongly broadened spectral lines in the optical) and the small radii of sdBs, the detection efficiency for companions via methods like radial velocity variations or transits is small. The timing of stellar pulsations offers a complementary detection method, sensitive to large orbital separations.
A small fraction of sdB stars shows pulsational variations in the p-(pressure-) and g-(gravity-) mode regimes. Rapid p-mode pulsators (sdBV r ), discovered by Kilkenny et al. (1997), show periods of the order of minutes and amplitudes of a few tens of mmag. Such pulsations were predicted by Charpinet et al. (1997, et seq.) to be driven by the κ-mechanism due to a Z-opacity bump. For slow pulsators (sdBV s ) the periods range from 30 to 80 min with small amplitudes of a few mmag. This class was discovered by Green et al. (2003) and the pulsations are explained by the κ-mechanism as well . Some sdB stars show both types of pulsation modes simultaneously (sdBV rs ). These hybrid pulsators lie at the temperature boundary near 28000 K between the two classes of pulsating stars, for example, the prototype for this class DW Lyn (Schuh et al. 2006), which is also addressed in this work, or Balloon 090100001 (Baran et al. 2005).
Pulsations driven by the κ-mechanism are coherent, which qualifies these objects for the timing method to search for substellar companions. This method is based on the light-travel time effect, with the host star acting as a stable 'clock' Spatial movements of the star around the barycentre induced by a companion result in time delays of the stellar light measured by the observer. Examples of detections using this method are 'pulsar planets' (e.g. Wolszczan & Frail 1992), planets detected by transit timing variations (e.g. Kepler 19 c, Ballard et al. 2011), planets orbiting δScuti stars (Murphy et al. 2016), or eclipsing binaries (e.g. V2051 Oph (AB) b, Qian et al. 2015). In particular, the detection of a late-type main sequence star companion to the sdB CS 1246 by Barlow et al. (2011b), subsequently confirmed with radial velocity data (Barlow et al. 2011a), or other studies like Otani et al. (2018), demonstrate the viability of this method in sdB systems. On the other hand, the particular example of V391 Peg b is currently under discussion (Silvotti et al. 2018) because of possible non-linear interactions between different pulsation modes that change arrival times (see Zong et al. 2018 for a detailed study of amplitude/frequency variations related to non-linear effects). Stochastically driven pulsations, are suspected by Reed et al. (2007a); Kilkenny (2010) and their nature confirmed by Østensen et al. (2014). Also, the candidate detections of KIC 05807616 and KIC 10001893 are uncertain, since other sdBs observed within the Kepler K2 mission exhibit g-modes with long periods up to a few hours. They question the interpretation of the low-frequency variations for KIC 05807616 and KIC 10001893 (e.g. Krzesinski 2015;Blokesz et al. 2019).
In order to detect sub-stellar companions orbiting rapidly pulsating sdB stars, the EXOTIME observational programme (EXOplanet search with the TIming MEthod) has been taking long-term data since 1999. EXOTIME conducted a long-term monitoring programme of five rapidly pulsating sdB stars. V391 Peg has been discussed by Silvotti et al. (2007Silvotti et al. ( , 2018. In this paper, we present the observations of DW Lyn and V1636 Ori, previously discussed in Lutz et al. (2008a);Schuh et al. (2010); Lutz (2011);Lutz et al. (2011), and re-evaluate their findings using an extended set of observations. In addition, the observations of QQ Vir and V541 Hya are presented and analysed. In the beginning of the programme, the mode stability was tested for all targets over a timespan of months in order to ensure the pulsation modes were coherent.
For the DW Lyn observations, Lutz et al. (2011) found no significant signals in a periodogram of the O − C data of the two analysed pulsation frequencies, which would indicate substellar companions. A tentative signal in the second frequency (in this work labelled f 2 , as well), formally corresponding to an 80-day companion orbit, is concluded to arise from mode beating of an unresolved frequency doublet. The analysis of V1636 Ori revealed a signal at 160 d in the periodogram of the main frequency O − C data (Lutz et al. 2011). Although this periodicity showed a significance of only 1 σ, Lutz et al. (2011) predicted an increase of significance with follow-up observations. We are using this extended data set in our work, now incorporating observations up to 2015.
This paper is organised as follows. Section 2 describes the observational aspects within the EXOTIME programme and the data reduction, followed by a description of our analysis in Section 3. Our results are presented in Section 4, together with a discussion.

Observations and data reduction
The observational data necessary for the analysis are comprised of many individual data sets gathered over the course of up to two decades. The detection method demands the observation of a target for a total time base at least as long as one orbit of a potential companion, which can span several years. This requires coordinated campaigns with observatories using~1 to 4 m telescopes. In order to derive sufficient accuracy for the analysis, observations with at least three to four consecutive nights, each with a minimum of two to three hours per target are required. To resolve the short-period p-modes the cadence must be shorter than about 30 s but still with a sufficient signal-to-noise ratio (S/N). All observations used the Johnson-Bessel B band. The correct time stamps for each observation are of most importance for the timing analysis. Most observatories of this study already successfully contributed to the work of Silvotti et al. (2007); Silvotti (2008, Table 2). The following list features some references where telescopes used for this study contributed successfully to other timing-relevant observations. Konkoly RCC 1.0 m Telescope: Provencal et al. (2009);Stello et al. (2006); Mt. Lemmon Optical Astronomy Observatory: Bischoff-Kim et al. References.

DW Lyn
There are photometric data available from 1999. Large gaps make a consistent O − C analysis difficult. Regular monitoring within the EXOTIME programme ran from 2007 until the beginning of 2010. Further observations cover a period up to the end of 2010. These multi-site observations are described in Lutz et al. (2008aLutz et al. ( ,b, 2011. Here, we add observations made with the SARA-KP 0.9 m telescope at Kitt Peak National Observatory in Arizona, that used exposure times of 30 s. Østensen et al. (2001) discovered V1636 Ori (HS 0444+0458) as a pulsating sdB star. Reed et al. (2007b) conducted a frequency analysis, reporting one small and two large amplitude p-modes.

V1636 Ori
V1636 Ori was observed between August 2008 and January 2015 for the EXOTIME project. About a third of the data was obtained using the 1 m South African Astronomical Observatory (SAAO) with the UCT and STE3 CCD instruments. Observations at the 1.2 m MONET/North telescope, equipped with an Apogee 1k×1k E2V CCD camera, were taken in 2x2 binnings, using 20 s exposure times. Observations at the 2.2 m Calar Alto Observatory (CAHA) used the CAFOS instrument with 10 s exposure time. Two nights were obtained at the 1.5 m telescope at Loiano observatory, using the BFOSC (Bologna Faint Object Spectrograph & Camera) instrument and 15 s exposure times. Between October 2008 and December 2009, observations at the 1 m Mt. Lemmon Optical Astronomy Observatory (LOAO) were conducted with a 2k×2k CCD camera with exposure times of 12 s and 20 s. The observations at the 3.6 m Telescopio Nazionale Galileo (TNG) in August 2008 and 2010 were performed with the DOLORES instrument and 5 s exposure times.
Observations of QQ Vir in 2001 and 2003 are described in Silvotti et al. (2002) and Silvotti et al. (2006), respectively. Between March 2008 and April 2010, the object was observed as part of the EXOTIME project (Benatti et al. 2010). Additionally, one observation run in February 2005 was performed at the 1.5 m telescope at Loiano observatory, using the BFOSC instrument.

TESS observations
The primary goal of the NASA Transiting Exoplanet Survey Satellite (TESS) space telescope is to detect exoplanets transiting bright nearby stars (Ricker et al. 2015). However, the extensive time series photometry is valuable for asteroseimology and the TESS Asteroseismic Science Consortium (TASC) coordinates short cadence observations of pulsating evolved stars. TESS observed V1636 Ori and V541 Hya with a cadence of 120 s between November 15, 2018 and December 11, 2018, and February 2, 2019 and February 27, 2019, respectively. We used the light curves provided by the MAST archive 1 that had common instrumental trends removed by the Pre-Search Data Conditioning Pipeline (PDC, Stumpe et al. 2012). Light curves and amplitude spectra are presented in Fig. A.1 and A.2. The twominute ('short') cadence undersamples the p-modes at about 140 s. In combination with the large photometric scatter, the amplitude spectra show no evidence of the p-modes. Thus, we did not make use of the TESS observations in our study.

Data reduction
For the EXOTIME observations, the data reduction was carried out using the IDL software TRIPP (Time Resolved Imaging Photometry Package, see Schuh et al. 2000). TRIPP performs bias-, dark-, flat-field corrections, and differential aperture photometry to calculate the relative flux of a target with respect to one  or more comparison stars and extinction corrections (second order polynomial in time). In the presence of sub-stellar companions, we might expect variations in the arrival times of stellar pulsations on the order of seconds to tens of seconds. The corresponding uncertainties are expected to be about one second. These uncertainties rise from observational constraints, such as smearing and sampling effects due to the integration time. The accuracy of individual time stamps is better than ±0.5 s. All time stamps were converted from GJD(UTC) to BJD(TDB), according to Eastman et al. (2010), with an accuracy well below the expected observational uncertainty. Typical S/N for our ground-based observations range from 60 for large amplitude pulsations, to 3 for the smallest pulsation amplitudes we investigate in this work. The amplitude spectra in Section 4 also show pulsations with smaller S/N, but these are not suitable for timing analysis because the uncertainties are too large (see Table B.1).

Analysis
In order to detect variations in the arrival time of stellar pulsations, we developed a pipeline to process the reduced data. A schematic flowchart of our pipeline is presented in Fig. 2. The input consists of the light curve (time series) and the dates of the observational epochs. In a light curve, typically spanning several years at a very low duty cycle of 0.2 to 1.7 per cent, an epoch consists of a few roughly consecutive nights of observation.
Outlier removal: In case no uncertainty in the flux measurement F is provided, the root-mean-square of each observation is used as an approximate photometric error for the later analysis. We have used a running median filter to exclude 5 σ flux-outliers. The length of the window size depends on the cadence of the observations. We constrained it to be not longer than half of the period of the main frequency. The analysis is performed for each frequency individually before all frequencies were analysed simultaneously.
Full data fit: For the individual fitting of pulsations, we first determined the frequency of the main signal. For this, we used the astropy package to calculate the Lomb-Scargle periodogram (Astropy Collaboration et al. 2013. From this periodogram, we selected the frequency with the largest amplitude to continue. In the next step, we performed the fit of a sinusoidal function to the light curve, using with amplitude A, frequency f , phase φ and offset o. The minimisation problem is solved using the scipy implementation of the Trust Region Reflective algorithm (Jones et al. 2011). The selected frequency from the amplitude spectrum serves as initial value, the full width at half maximum of the corresponding peak in the periodogram is used as a boundary. The amplitude-guess is taken from the amplitude spectrum. In case of highly varying amplitudes, the initial value can be set manually. The fitting routine returns the parameters and their variance.  Epoch fit: Frequency and offset are all kept fixed for the following analysis of the individual observational epochs. The starting value of the current phase fit is determined by the average of the previous j phase values (or the global fit value from above in case there are no j previous values yet) in order to keep the fitting process stable and avoid 'phase-jumps'. For our target sample, a value of j = 3 has proven to be reasonable, except when observational gaps span over several years. The uncertainty in the phase measurement scales inverse with the length of the epochs. Thus, this length is chosen in a way to minimize the uncertainties of the fit but at the same time keep the epochs as short as possible to maximize the temporal resolution of the final O − C diagram. Often, the observations themselves constrain the length of the epochs (e.g. three consecutive nights of observations and a gap of several weeks before the next block of observations). If possible, we aimed for an epoch length such that the timing uncertainties are of the order of one second. The phase information of the global and the epoch fit As a last step in the single-frequency analysis, the fitted model is subtracted from the light curve. We noticed significant amplitude variations for some of our targets. Thus, we subtracted the model using the amplitude of the individual epochs. This prewhitening procedure is repeated for every relevant pulsation in the data.
Multi frequency fit: Close frequencies are likely to introduce artificial trends in the arrival times in such a step-by-step analysis. Thus, the sum of all sinusoidal functions, is fitted to the non-whitened light curve, where n is the number of investigated frequencies. The previously retrieved values for amplitude, frequency, phase and offset are used as initial values. We used the phase information φ n of the light curve as reference phase, namely calculated phase C in the final O − C diagram. Similar to the single-frequency analysis, the observational epochs are fitted individually using the sum of sinusoidal functions to yield the observed phase information O.
The results of the simultaneous fit for each target in this paper are summarised in Table 3. We list pulsation modes not used for the timing analysis in Table B.1. Fig. 3, 4, 5, and 6 show example light curves of the targets for one epoch each, including their multi frequency fit and the respective amplitude spectrum.

Results and discussion
In the following, we discuss the implications of the obtained amplitude spectra and O−C measurements on the evolutionary state and presence of sub-stellar companions to the targets. Table 3: Parameters of the simultaneously fitted pulsations per target over their full observational time span and the pulsation period P. The phase φ refers to the time corresponding to the first zero-crossing of the function after the first measurement t 0 in MBJD. 6.72 6.74 6.76 6.78 6.80 6.82 6.84 6.86 6.88 6.90

DW Lyn
The amplitude spectrum of DW Lyn in Fig. B.1 reveals two strong pulsation modes at f 1 = 237.941160 d −1 and f 2 = 225.15898 d −1 . A closer look to the amplitude spectrum in Fig. 7 reveals small asymmetries compared to the window function.
The pre-whitening of both frequencies leaves residuals well above noise level in the amplitude spectrum, indicating unresolved multiplets or mode splitting, especially for f 2 . The S/N of modes at higher frequencies, for example, at about 320 d −1 and 480 d −1 , are too small for a stable O − C analysis (see Table B.1). Therefore, the O − C diagram in Fig. 8 shows the analysis of the two main pulsation modes, with the time-dependent variation of the pulsation amplitudes.
In order to determine evolutionary timescales of the pulsations, we investigated the long-term evolution in the O − C data. A constant change in period results in a second-order term as a function of time (Sterken 2005), which allows us to derive a value for the secular change of the periodṖ, and hence the evolutionary timescale. Results of the fits of the second order polynomial are included in Fig. 8, which areṖ/P f 1 = (5.8 ± 0.2) × 10 −5 d −1 andṖ/P f 2 = (−29.3 ± 0.8) × 10 −5 d −1 . AssumingṖ is based on stellar evolution, stellar model calculations show that the sign of the rate of period change indicates the phase of the sdB after the zero-age extreme horizontal branch (ZAEHB) (Charpinet et al. 2002). For p-modes, a positiveṖ relates to the first evolutionary phase of the ZAEHB, in which the surface gravity decreases due to He burning in the core. A neg-ativeṖ would correspond to the second evolutionary phase, in which the sdB contracts because the depletion of He in its core, and this happens before the post-EHB evolution. The turning point between these two states occurs between 87 and 91 Myr after the ZAEHB. According to our measurement of a positiveṖ for f 1 , DW Lyn would still be in its first evolutionary phase. With the lack of a mode identification from an asteroseismic model for DW Lyn, we can not directly compare the measuredṖ with theoretical predictions from Charpinet et al. (2002)  P with a comparable order of magnitude to our measuremenṫ P = (4.3 ± 0.15) × 10 −1 s Myr −1 , for example,Ṗ = 1.62s Myr −1 for a model with a mode of l = 0, k = 0 at the age of 67.83 Myr (Charpinet et al. 2002, appendix C). The largeṖ of f 2 is consistent with the apparent mode splitting seen in the amplitude spectra in Fig. 7, and thus does not reflect the evolutionary phase of DW Lyn. After subtracting the long-term trend, small timescale features are evident. For example, the O − C data for f 2 show an oscillating behaviour with a significance of 3 σ within the first 200 days, while the arrival times for f 1 remain constant during the same period of time. In later epochs, the O − C data for both frequencies agree mostly within 2 σ. During the second half of the observations, the phase of f 2 jumps by about 100 s. This behavior lacks an explanation.
Additionally, the evolution of the pulsation-amplitudes in Fig. 8 shows a comparable oscillating behaviour for f 2 within the first epochs similar to the change in arrival times. Although the periodic variations in amplitude are not as significant as for the phase, the occurrence of simultaneous phase-and amplitudemodulations indicate a mode beating of two close, unresolved frequencies. The residuals in the amplitude spectrum support this explanation. In later observations, the amplitude remains almost constant within the uncertainties. The beating mode might lose energy or shift frequency over time. The amplitude for the f 1 pulsation drops by about 1 per cent (amplitude), or about 35 per cent (relative) to the second half of the observation campaign with a similar quasi-periodic variation as the phase. The residuals in the amplitude spectrum show no indication of an unresolved frequency leading to mode-beating. Besides stochastically driven pulsation modes, Kilkenny (2010) suggested energy transfer between modes as possible explanation for amplitude variations. For both frequencies, a possible interaction between amplitude and phase of pulsations is not well understood.

V1636 Ori
The amplitude spectrum of V1636 Ori in Fig. B.2 shows two main pulsation modes with frequencies at f 1 = 631.7346 d −1 and f 2 = 509.9780 d −1 . The S/N is not sufficient to use a third pulsation mode at 566.2 d −1 (6553 µHz, Reed et al. 2007b). The amplitude spectrum of TESS data in Fig. A.2 shows no evidence for g-mode pulsations with amplitudes greater than 0.4 per cent. A detailed look at the spectra of the two main frequencies in Fig. 9 shows mode splitting, likely due to a change in frequency over the long observation time.
The O − C diagram in Fig. 10 shows the two main pulsation modes and the variation of the pulsation amplitudes.
From the second order fit in time, we derive the changes in periodṖ/P f 1 = (−8.54 ± 0.14) × 10 −5 d −1 andṖ/P f 2 = (−2.5 ± 0.5) × 10 −5 d −1 . We caution the interpretation of these values as evolutionary timescales since the apparent mode splitting seen in Fig. 9 could explain these trends as well.
The residuals after subtracting the long term trend show a large variation. They change by up to about ±50 s for f 1 (∼ 14 σ significance) and up to about ±30 s for f 2 (∼ 3 σ significance). The amplitude for f 1 drops by about 0.25 per cent (amplitude) or about 33 per cent (relative) in the time between MBJD = 55100 d and 55300 d, and returns to its previous level afterwards, while the amplitude for f 2 remains constant within the uncertainties. This decrease in amplitude coincides with earlier arrival times in the O − C diagram. As already discussed in the previous section, a possible amplitude-and phase-interaction is not well understood. The f 1 pulsation mode may not be coherent on such long timescales but of a short-term stochastic nature not resolvable by our data set (e.g. KIC 2991276, Østensen et al. 2014).  Fig. 11 in detail and shows asymmetries compared to the window function. After the pre-whitening process, a close frequency at about 626.881270 d −1 remains but attempts to model this pulsation fail with uncertainties too large for the timing analysis. There appear two more frequencies suitable for our study. The amplitude spectra around f 2 = 552.00713 d −1 and f 3 = 642.0516 d −1 are presented next to f 1 in Fig. 11. An-other peak at about 665 d −1 consists of at least two frequencies at 664.488549 d −1 and 665.478133 d −1 , but they are not sufficiently resolvable within the individual epochs, and lead to uncertainties in the O − C analysis that are too large. Figure 12 shows the resulting O − C diagram and the amplitudes at different epochs. Due to the large observational gap from 2003 to 2008 with only one block of observations in between, we had difficulties avoiding errors in cycle count. In order to avoid a phase jump, we increased the averaging window for initial phase values to q = 6. With this set up, the changes in pulsation frequencies read as follows:Ṗ/P f 1 = (1.7 ± 1.6) × 10 −7 d −1 , Article number, page 8 of 22 P/P f 2 = (2.4 ± 0.4) × 10 −5 d −1 andṖ/P f 3 = (4.0 ± 0.5) × 10 −6 d −1 . While f 2 and f 3 show no significant variation of pulsation amplitude, f 1 varies by 1.5 per cent (amplitude) or 50 per cent (relative). Thus, the corresponding phase changes should be interpreted with caution. Charpinet et al. (2006) identified the radial order k and degree l from asteroseismic modelling to be f 1 : l = 2, k = 2; f 2 : l = 4, k = 1; f 3 : l = 3, k = 2. These combinations do not allow a direct comparison of ourṖ measurements to the model calculations from Charpinet et al. (2002), but the sign ofṖ indicates QQ Vir to be in the stage of He burning.

V541 Hya
The amplitude spectrum in Fig. B.4 shows two pulsation modes with frequencies at f 1 = 635.32218 d −1 and at f 2 = 571.28556 d −1 . Both of them show a complex behaviour (Fig. 13), indicating unresolved multiplets and/or frequency changes that we see also in the O-C diagrams (Fig. 15). The S/N for a third frequency at 603.88741 d −1 is not sufficient for the O − C analysis. Similar to V1636 Ori, the amplitude spectrum obtained from the TESS light curve in Fig. A.2 shows no evidence for g-mode pulsations with amplitudes greater than 0.4 per cent. Randall et al. (2009) speculated about rotational mode splitting for f 3 with ∆ f 3,− = 5.12 µHz and ∆ f 3,+ = 3.68 µHz. The asteroseismic modelling associates f 1 with a l = 0 mode and f 2 with l = 0 or 1 mode (depending on the favoured model). f 3 corresponds to a l = 2 mode. They caution this interpretation due to their limited resolution in frequency space, the mode splitting could be an unresolved quintuplet. Our data set shows no clear evidence for a mode splitting with ∆ f 3,− = 5.12 µHz or ∆ f 3,+ = 3.68 µHz (see Fig. 14) but rather a mode splitting for f 1 and f 2 with about ∆ f = 0.08 µHz (Fig. 13). Assuming these modes are of degree l = 1, this could be interpreted as a triplet. But Randall et al. (2009) model these modes with a degree of l = 0, which does not support a mode splitting into triplets.
The O − C diagram in Fig. 15 shows the analysis of the two main pulsation modes and the variation of the pulsation amplitudes. The second order fits in time correspond to changes in period ofṖ/P f 1 = (−1.49 ± 0.11) × 10 −5 d −1 andṖ/P f 2 = (−0.7 ± 1.5) × 10 −5 d −1 . For f 2 , the change in period does not significantly differ from the null hypothesis. Assuming these changes origin from stellar evolution, V541 Hya might just have passed the point of sign change inṖ and at the beginning of the contraction phase. While the arrival times scatter widely, the amplitudes of both pulsations remain almost constant within the uncertainties. If V541 Hya is in its evolution close to starting the contraction phase, as indicated by aṖ close to zero, the changes in stellar structure may cancel the strict phase coherence.

Testing the sub-stellar companion hypothesis
In order to set upper limits to the mass of a companion, we computed a series of synthetic O − C curves for different orbital periods and companion masses, assuming circular orbits, and compared these curves with the O−C measurements after subtracting the long-term variations. For each synthetic O − C curve, we selected the phase that gives the best fit to the data using a weighted least squares algorithm. For each observational point, we computed the difference, in absolute value and in σ units (where σ is the O − C error), between O − C and the synthetic value. The greyscale in Figs. 16 and 17 corresponds to the mean value of this difference in σ units, which means that the presence of a companion is indicated by a minimum (bright areas) of this parameter. We see that in V1636 Ori, QQ Vir, and V541 Hya, the mean difference for f 1 is always very high, implying that the data are not compatible with a companion. However, these results are limited by the fact that the O −C diagrams of these stars are 'contaminated' by other irregular variations, presumably due to other reasons like non-linear interactions between different pulsation modes, for example, and therefore these constraints to the orbital period and mass of a companion must be taken with some caution. For the f 2 and f 3 measurements, the mean difference to the synthetic data is smaller in sigma units (because of the larger uncertainties) and very uniform. The uncertainties of the O-C measurements are not small enough to favour a set of models in the period-mass parameter space.
For f 1 of DW Lyn, there is a significant minimum at about 1450 d (∼ 4 years) and ∼ 5 M 2 , which is also well visible in the O − C diagram of Fig. 8. This periodicity is not visible in the second frequency f 2 which, however, has much larger error bars due to the much lower amplitude of f 2 with respect to f 1 . Lutz et al. (2011) described a periodicity at 80 days, detected for f 2 . We can recover this signal, however, with a low significance. This would correspond to a light-travel time amplitude of 4 s (for m sin i ≈ 15 M ), which is smaller than the 15 s measured by Lutz et al. (2011). Nevertheless, this signal is not confirmed by f 1 . Thus, we rule out a companion induced signal in the arrival times due to the lack of simultaneous signals in f 1 and f 2 with similar amplitude. The tentative signal in f 2 is better explained by mode beating, as already described in Sect. 4.1. The variations seen in the first 200 days of the O − C diagram in Fig. 8 correspond to a periodicity of about 80 days and are accompanied by variations in the amplitude of the pulsation.
For V1636 Ori, Lutz (2011) predicted a period at 160 d and amplitude of 12 s. This can not be confirmed as a companioninduced signal. A periodic signal with an amplitude of 6.5 s (for m sin i ≈ 15 M ) is indicated in the analysis of f 1 , but at a low significance and accompanied by many other signals of similar significance. This periodicity is not confirmed by a significant signal in the measurements of f 2 .

Summary and conclusion
In this work, we present ground-based multi-site observations for the four sdBs, DW Lyn, V1636 Ori, QQ Vir, and V541 Hya. We investigated variations in the arrival times of their dominant stellar pulsation modes to draw conclusions about secular period drifts and possible sub-stellar companions. All light curves are analysed homogeneously. From the O − C measurements, we derive an evolutionary timescale from the change in periodṖ. Comparing to model calculations from Charpinet et al. (2002), we infer the evolutionary phase of the target. Although someṖ measurements are influenced by mode splitting, we can tell from the sign ofṖ 1 of DW Lyn that the star is likely still in the stage of central He burning. We can draw a similar conclusion from the sign ofṖ of QQ Vir. TheṖ measurements of V1636 Ori are likely affected by mode splitting, making it difficult to interpret the results in the context of stellar evolution. V541 Hya showsṖ measurements close to zero, which indicates the star being at the transition phase between He burning and contraction due to the depletion of He in the core.
Comparing the atmospheric properties from Table 1 with the evolutionary tracks for different models from Fig. 1 in Charpinet et al. (2002), we can confirm the hypothesis that DW Lyn and QQ Vir are in their He burning phase. V541 Hya agrees within 2 σ of the log g measurement with one model at the turning point between the two evolutionary stages.
However, we can not exclude frequency and amplitude variations on smaller timescales than resolvable by our data set. Using temporally higher resolved Kepler-data of KIC 3527751, Zong et al. (2018) cautioned about long-term frequency or phase evolutions ascribing to non-linear amplitude and frequency modula-tions in pulsating sdBs. We see such effects already in our data set, even with a low temporal resolution compared to the Kepler sampling with a duty cycle of more than 90 per cent.
Observations on DW Lyn and V1636 Ori were published by Lutz et al. (2008a);Schuh et al. (2010); Lutz (2011);Lutz et al. (2011). Our analysis of these observations, including extended data sets, do not confirm the tentative companion periods of 80 and 160 days, respectively. These signals more likely arise due to mode beating indicated by partly unresolved frequency multiplets and amplitude modulations.
Almost all analysed pulsation modes show formal significant changes in arrival times, but the amplitudes of these periodic signals do not correlate with frequencies, excluding the lighttravel time effect due to orbital reflex motions for such variations and thus giving upper limits on companion masses. Only DW Lyn might have a planetary companion on a long orbital period, as indicated by one arrival time measurement. But this can not be confirmed with a second measurement, due to larger uncertainties. Additionally, more studies question the presence of already proposed companions, for example, Krzesinski (2015); Hutchens et al. (2017). Our unique sample of long-term observations shows a complex behaviour of mode-and amplitude interactions in sdBs which should be addressed in further studies. Until this has been addressed, caution is advised when interpreting O − C pulse arrival times in terms of companions.       (6) 0.01(7) 502.410 (2) 0.01(9) V541 Hya 531.16759(16) 0.03 (7) 603.88741 (6)