A new photometric and dynamical study of the eclipsing binary star HW Virginis Free

ABSTRACT

A growing number of eclipsing binary systems of the ‘HW Virginis’ (HW Vir) kind (i.e. composed by a subdwarf-B/O primary star and an M dwarf secondary) show variations in their orbital period, also called eclipse time variations (ETVs). Their physical origin is not yet known with certainty: While some ETVs have been claimed to arise from dynamical perturbations due to the presence of circumbinary planetary companions, other authors suggest that the Applegate effect or other unknown stellar mechanisms could be responsible for them. In this work, we present 28 unpublished high-precision light curves of one of the most controversial of these systems, the prototype HW Vir. We homogeneously analysed the new eclipse timings together with historical data obtained between 1983 and 2012, demonstrating that the planetary models previously claimed do not fit the new photometric data, besides being dynamically unstable. In an effort to find a new model able to fit all the available data, we developed a new approach based on a global-search genetic algorithm and eventually found two new distinct families of solutions that fit the observed timings very well, yet dynamically unstable at the 10⁵-yr time-scale. This serves as a cautionary tale on the existence of formal solutions that apparently explain ETVs but are not physically meaningful, and on the need of carefully testing their stability. On the other hand, our data confirm the presence of an ETV on HW Vir that known stellar mechanisms are unable to explain, pushing towards further observing and modeling efforts.

1 INTRODUCTION

The discovery of the first exoplanets by Wolszczan & Frail (1992) and Mayor & Queloz (1995) was the starting point to the detection of a great number of other planetary systems through different observing techniques. Although the majority of them have been found orbiting Sun-like stars (e.g. Petigura, Marcy & Howard 2015), there is an increasing number of exoplanets being discovered orbiting all kinds of stars (e.g. Gould et al. 2014; Gillon et al. 2017; Brewer et al. 2018). A particularly interesting case among them is represented by circumbinary planets, which orbit a binary system instead of a single star. These kinds of planets can be detected, among other techniques (such as transits, e.g. Kostov et al. 2016; radial velocity, e.g. Konacki et al. 2009; or light traveltime delay, e.g. Silvotti et al. 2018), by measuring and analysing changes in the orbital period of eclipsing binary stars, a dynamical method commonly known as eclipse time variations (ETVs; e.g. Sale et al. 2020). These variations have been observed in a wide range of binary systems, such as post-common envelope binaries, which exhibit modulation periods of a few tens of years (e.g. Bours et al. 2016). A possible mechanism to explain ETVs is the light traveltime effect (LTTE; also known as the Rømer effect), which refers to the combination of the motion of the stellar components with respect to the barycentre of the system due to the gravitational perturbation of additional bodies, with the finite speed of light (Irwin 1952).

Among the vast taxonomy of eclipsing binaries, the so-called ‘HW Virginis’ (HW Vir) systems have recently drawn the attention of astronomers. These systems are post-common envelope binaries composed of a subdwarf of spectral type O or B and a late-type main-sequence star, e.g. sdB + dM for the prototype. They have very short orbital periods (of the order of a few hours), and in a surprisingly high fraction of the cases, ETVs have been observed, typically from tens of seconds to several minutes of amplitude and semiregular modulations on long time-scales, from years to decades (see Heber 2016, for a detailed review on HW Vir systems). Different explanations have been proposed to interpret ETVs, usually based on two different effects or a combination of them: the LTTE effect caused by one or more unseen companions, and the so-called Applegate effect. The latter was first proposed by Applegate (1992), and it interprets the variations on the orbital period as a consequence of magnetic activity in one of the stars of the binary system (in the case of HW Vir, the main-sequence component). According to Applegate (1992), the distribution of the angular momentum in the active star changes as the star goes through its activity cycle. These variations on the angular momentum distribution induce a change in the gravitational quadrupole moment of the star (making it more or less oblate), which can cause perturbations in the orbit of the system and thus in the orbital period.

In this work, we analyse data from the prototypical HW Vir, a detached eclipsing binary system first identified as such by Menzies & Marang (1986). HW Vir has a very short period of 2.8 h, and its components have masses of 0.49 and 0.14 M_⊙, for the sdB and dM components, respectively (see Table 1 for the most recent parameters of HW Vir). Since its discovery, the system has been broadly studied due to its intrinsic characteristics and its striking period variations. A decrease in the orbital period of the system was first detected by Kilkenny, Marang & Menzies (1994), followed by Çakirli  & Devlen (1999), who re-analysed the eclipse timings between 1984 and 1999 and concluded that LTTE was the most promising explanation for the observed period variations. They proposed that HW Vir was revolving about a third body with a period of 19 yr. Later on, further studies were performed (Wood & Saffer 1999; Kilkenny et al. 2000; Kiss et al. 2000) analysing the period variations with different techniques, without reaching a definitive explanation. Kilkenny, van Wyk & Marang (2003) presented new eclipse timings for HW Vir and confirmed the presence of a periodic LTTE term due to a third body in the system, a claim also supported by İbanoǧlu et al. (2004).

Lee et al. (2009) presented new CCD photometry with an 8-yr baseline, and proposed that the linear term of the period decrease (dP/dt) may be caused by angular momentum loss due to magnetic stellar wind braking, while the cyclic period variations may be interpreted as LTTE terms induced by the presence of two additional bodies in the system, having masses of M₃sin i₃ = 19.2 M_J and M₄sin i₄ = 8.5 M_J, respectively.¹ This model was independently tested by Beuermann et al. (2012), who found that it fails to fit their new eclipse timings and it is dynamically unstable on a time-scale of a few thousand years. Beuermann et al. also proposed a new LTTE model with two companions with masses of M₃sin i₃ ≃ 14 M_J and M₄sin i₄ = 30-120 M_J, and periods of 12.7 and 55 ± 15 yr, respectively. Horner et al. (2012) independently tested Lee et al.’s model and came to the same conclusion about the dynamical instability of the system on very short time-scales; they also claimed that the ETVs cannot be driven by gravitational influence of perturbing planets only, and that there must be another astrophysical mechanism taking place in order to explain them.

Finally, Esmer et al. (2021) found a new two-planet solution, but it did not appear to be dynamically stable. The main differences between our approach and theirs will be summarized in the ‘Discussion and Conclusions’ section.

Regarding the Applegate effect, Navarrete et al. (2018) analysed the required energy to drive the Applegate effect in a sample of 12 close binary systems (including HW Vir), and compared it with the energy production of a simulated sample of magnetically active stars. In the case of HW Vir, they discarded the possibility of this effect being the underlying cause for the ETVs, since the magnetic field of the magnetically active star (i.e. the dM star) is not strong enough to produce these variations.

A conclusive explanation to HW Vir’s ETVs is still missing. For this reason, our aim is to derive new eclipse timings from our unpublished photometric data, and use them along with the ones available in the literature to better constrain the physical parameters characterizing the system of HW Vir, as well as to test these new parameters for dynamical stability on a large time-scale.

The paper is organized as follows: In Section 2, we present our data, along with the data reduction process we followed, the light-curve fitting, and the determination of the eclipse timings, while in Section 4 we outline the LTTE modeling and test the previous model proposed to explain the ETVs of the system with our new data, as well as using an N-body integrator to test its dynamical stability. In Section 5, we describe the method we used to estimate new parameters for the putative companions of HW Vir. In Section 6, we discuss our findings and we draw some conclusions regarding the explanation behind the ETVs of HW Vir as well as some prospects for future work.

2 OBSERVATIONS AND DATA REDUCTION

Our analysed data set consists of 30 photometric observations of HW Vir obtained over time span of ∼11 yr (2008–2019), including 28 previously unpublished light curves. For our analysis, we combined data from five different instruments as described in the following.

From the Asiago Astrophysical Observatory located on Mt. Ekar in Asiago, we obtained 15 light curves using the |$1.82\,\rm{-m}$| ‘Copernico’ telescope and the Asiago Faint Object Spectrograph and Camera. These images were taken with an exposure time ranging from 2 to 6 s, through the V, R, and r filters. Three light curves were obtained using the |$67/92\,-\,\rm{cm}$| Schmidt telescope located at the same observatory. These observations were carried out in the R and r filters with an exposure time of 20 s, except for the last one (4 s).

Six light curves were obtained using the telescopes of the ‘Gruppo Astrofili Salese Galileo Galilei’;² the telescopes have a primary mirror with a diameter of 410 mm and a focal length of 1710 mm and they are located in Santa Maria di Sala, in Northern Italy. The observations were carried out in the V filter with an exposure time ranging from 20 to 45 s.

Our largest set in terms of number of data points comes from the Wide Angle Search for Planets (WASP-South; Pollacco et al. 2006), a transit survey with an array of small telescopes operating at the South African Astronomical Observatory (SAAO). WASP-South gathered four full seasons of observations of HW Vir from 2008 to 2012, for a grand total of 353 measured primary eclipses. This particular data set has not yet been included in a public data release, and has been kindly provided to us by the WASP-South team.

We also include in our analysis two light curves from K2 (Howell et al. 2014), observed during Campaign 10, and a vast collection of literature timings already analysed by Beuermann et al. (2012) and summarized at the end of this section. A detailed summary of all the observations is given in Table 2. Each light curve is identified with a unique ID with the leading letter matching the telescope: w for WASP-South, s for Asiago Schmidt, g for GAS, c for Asiago Copernico, and kt for K2. The w and kt light curves are split into four and two ‘chunks’ (respectively), for the reasons explained in Section 3.

Due to the lack of stellar crowding in the field of HW Vir, we use the differential aperture photometry technique to reduce our photometric series from the c, s, and g data sets. To perform the usual data reduction and the aperture photometry, we use the software STARSKY, a pipeline written in FORTRAN 77/90 by Nascimbeni et al. (2011, 2013), that was specially developed for The Asiago Search for Transit timing variations of Exoplanets (TASTE) project. As for the w data set, we take the light curves as they were delivered by the standard WASP software pipeline. For the K2 data, we extracted the light curve by reconstructing the 89 970 images containing HW Vir as done in Libralato et al. (2016), and performing a three-pixel aperture photometry of the target on each image, subtracting the local background measured in an annulus centred on the target and having radii r_in = 7 pixels and r_out = 15 pixels to the total flux. We detrended the light curve following the procedure by Nardiello et al. (2016). The resulting light curves from all the observations are shown in Fig. 1.

In order to measure timing variations with an absolute accuracy much better than 1 min, as needed for measuring ETVs, it is crucial to convert all our timestamps to a single, uniform time standard. Therefore, we convert all of them to the so-called Barycentric Julian Date computed from the Barycentric Dynamical Time, or BJD_TDB, following the prescription by Eastman, Siverd & Gaudi (2010). For this task, we rely on the VARTOOLS code.³ Due to the crucial importance of this step for our dynamical analysis, we perform a double check of the conversion with the help of the online tool⁴ made available by Eastman et al. (2010). We also apply this time conversion to all the 287 literature timings from SAAO, Wood, Zhang & Robinson (1993), Lee et al. (2009), BAV, VSNET, AAVSO, BRNO, and Beuermann et al. (2012), who, in turn, used timings from MONET/North. Again, all the HJD_UTC and BJD_UTC are homogeneously converted to BJD_TDB to ensure a proper comparison between the old timings and our new ones. A comprehensive listing of all the literature timings as converted by us is available in Appendix  A.

3 ECLIPSE TIMINGS

To retrieve the best estimate of the orbital and physical parameters of the system, and most crucially the eclipse central time T₀, we fit an appropriate model to our light curves. For this purpose, we use the JKTEBOP ⁵ code (Southworth 2012), which was originally developed to fit light curves of detached eclipsing binaries and later adapted to model also exoplanetary transits. JKTEBOP implements non-linear least-squares optimization techniques [based on the Levenberg–Marquardt (LM) algorithm; Moré 1978]. It has different ‘tasks’ to choose from, according to how the light curves would be fitted and how the uncertainties are estimated. This process is meant to determine the best-fitting values of T₀ for each individual light curve and a reliable error estimate.

As a first step, we check that the software is properly fitting our light curves and converging to a physical solution by using task3, i.e. by simply running the task to each preliminary light curve and performing a visual inspection. At this stage, we decide to split the w and kt light curves into separate ‘chunks’. For the WASP-South data, this is done because the composite light curve has a 4-yr coverage, and fitting it as a whole could in principle smear the LTTE signal; by splitting it into four distinct ‘seasons’ of about 4 months each, we completely avoid this risk (the shortest significant O − C periodicity reported in the literature being ∼3000 d). As for the K2 data, the Campaign 10 light curve shows a large 2-week gap due to a repointing procedure followed by an unexpected shutdown of the camera. To make ourselves sure that there are no systematic errors introduced by this issue, we separately analysed the two chunks before and after the blank gap.

We then remove the outliers from our light curves at 4σ using task4 of JKTEBOP, and, since we want to obtain a reliable measure of the eclipse time (T₀), we need to first build consistent templates of the parameters for each of the filters of our observations, to leave only T₀ as a free parameter in the final fit. To do this, we join the full-phase light curves from the same filter (since the light curves are colour dependent) and leave the following parameters free to find the best-fitting values: the sum of the stellar radii (R₁ + R₂), their ratio R₁/R₂, the inclination of their orbit, the surface brightness and the limb darkening of the primary star, the reflection coefficient of the secondary star, the scale factor, and the eclipse time (T₀). We do this for the V and R/r filters, and additionally, for the WASP and K2 light curves. Then, we run task9 of JKTEBOP, which uses a residual-shift method to obtain the best fit. This method evaluates the best fit for the data points and shifts the residuals of the fit point by point through all the data, calculating a new best fit after each shift. This approach allows us to have as many best fits as points in the input light curve, and it also estimates the relevance of the correlated red noise to the parameters of the fit. The output of this task is therefore three high-accuracy parameter sets (templates), one for each filter: a V template for the Copernico/V and GAS light curves; an R/r template for the Copernico/R, r and the Schmidt/R light curves; and an unfiltered template for the WASP and K2 light curves.

We retrieve the T₀s by running task9 one more time, fixing all the parameters except the eclipse times. An example of the quality of the fit on our two most complete light curves from the Copernico telescope (⁠|$\tt {c2}$| in Bessel V and |$\tt {c11}$| in Sloan r) is shown in Fig. 2.

The resulting timings of HW Vir are reported in Table 3. We compute a total of 30 mid-eclipse timings, with an excellent median timing error for our light curves of only ∼1.3 s and down to 0.3 s for the best ones (from the c and w sets). Our new data increase the current number of high-precision observations [σ(T₀) < 5  s] by about 50 per cent, and extend the baseline by 6 yr with respect to the dynamical study of HW Vir (Beuermann et al. 2012).

4 MODELLING

4.1 LTTE calculation

The approach of Irwin (1952) was to use the plane perpendicular to the line of sight that passes through the centre of the elliptical orbit of the binary about the centre of mass of all the bodies in the system as the reference frame, which adds a ne_ksin ω_k term to equation (2). Our approach is to use another perpendicular (and parallel) plane to the line of sight that passes through the centre of mass of all the bodies in the system as the reference frame, resulting in the exclusion of this term.

4.2 Test of the previous model

To test this two-companion model, we plot the O − C diagram using both the literature data and our new eclipse timings in Fig. 4. The model is able to reproduce the data from the literature very well; however, it fails to fit our new data.

We check the dynamical stability of this model by reproducing the same test performed by Beuermann et al. (2012) using the MERCURY6 ⁶ (Chambers 1999) package. We set the initial Keplerian parameters of the system with the binary as a single body of mass M_bin = M₁ + M₂ at the centre of the system, as described in Beuermann et al. (2012), and we use the same hybrid symplectic integrator. As a first test, we integrate for 10⁴ yr, and we find that the inner planet is ejected after ∼2500 yr, in contrast with Beuermann et al. (2012)’s paper, who suggest that their proposed model is stable for 10⁸ yr.

We perform additional checks using the radau integrator within the MERCURY6 code, and also using the PYTHON-C package rebound ⁷ (Rein & Liu 2012) with three of their different integrators namely, ias15 (Rein & Spiegel 2015), whfast (Rein & Tamayo 2015), and mercurius. All the simulations were run for 10⁶ yr, using a step size of 8.8 d (1/530 of P₃) with output every 308.9 d (1/15 of P₃). Additionally, we test the stability with a new version of MERCURY6, MERCURY6_binary,⁸ a modified version of the original code by Smullen, Kratter & Shannon (2016), which allows us to simulate both single and binary stars, treating the central star in the binary as a composite ‘big body’ instead of a single central object. Following the advise by the author, we use the radau integrator to perform the simulation, and we integrate for 10⁶ yr with the same step size described above. We consider a planet to escape or be ejected at a distance >150 au.

The initial orbital and physical parameters used for all the simulations performed are listed in Table 4. The results of all the simulations returned unstable systems, in different time-scales and for different reasons, such as ejection of outer or inner planet, a close encounter between planets, or the inner planet colliding with the binary. As a final check, we use the Mean Exponential Growth factor of Nearby Orbits (MEGNO; Cincotta & Simó 2000) indicator in rebound. Briefly, the MEGNO indicator 〈Y〉 will reach the value of 〈Y〉 = 2 for stable orbits, and it will be 〈Y〉 ≫ 2 for unstable configurations (in the case of 〈Y〉 > 4 or a close encounter and an ejection, we assign the maximum value 〈Y〉 = 4). We set the initial conditions as in Table 4, but we let vary, for the inner companion (identified with the subindex 3), the semimajor axis a₃ from 1 to 6 au and the eccentricity e₃ from 0 to 0.5, both in 100 linear steps. We compute the orbits of each configuration with the whfast integrator with a step size of 1 d for an integration time of 10⁵ yr. The final grid has 100 000 simulations, each returning a MEGNO value. As shown in Fig. 5, we find that the solution from Beuermann et al. (2012), depicted by the red dot, lies on an unstable region, confirming our tests with different codes and integrators. It is worth noting that all the simulations have the same reference frame as in Winn (2010), which is the plane X–Y in the sky plane and Ω_{3, 4} = 180°, and we assume the orbits to be coplanar with the binary.

5 A NEW MODEL

Our aim at this stage is to find a new LTTE model that properly fits the data. We separately analysed two data sets: one with all the available data (317 points), and one for which we discarded the first two observing seasons from the literature (35 photoelectric measurements between JD 2445730 and 2445745 from Kilkenny et al. 1994). From now on, we will refer to these data sets as the ‘full’ and the ‘reduced’ one, respectively. The latter selection was done as a test since the Kilkenny et al. (1994) data were always suspiciously offset from any best-fitting model and lack the original time-series data; i.e. we are unable to perform any independent check on them. We also rescale all the T₀ errors by adding in quadrature 1 s to Beuermann et al. (2012)’s and our values, and 5 s to the rest of the literature values. We apply this rescaling to take into account systematic errors in the absolute calibration of the timestamps at this level (due for instance to clock drift, to the finite shutter traveltime, or to technical dead times while commanding the camera or saving the images). This assumption will be later empirically justified by the residual of our best-fitting models being very close to |$\chi ^2_r\simeq 1$|⁠.

After removing the outliers and rescaling the errors, we extend the code described in Section 4.1 with the implementation of PIKAIA (Charbonneau 1995), a genetic algorithm to solve multimodal optimization problems. This algorithm is based on the theory of evolution by means of natural selection; that is, a new population is generated by choosing the fittest pairs from the original population, and this process continues until a certain fitness level is achieved or after a pre-defined number of generations. We perform 100 000 simulations of 1000 generations each on a population of 200 individuals and we use the inverse of the reduced chi-square |$1/\chi _\mathrm{r}^2$| as our fitness function. Once the code computes the results for PIKAIA at the end of each simulation, it uses the LM algorithm to refine the PIKAIA output and it calculates the final best-fitting solution.

We also run an independent analysis based on a modified version of PIKAIA in FORTRAN 90, wrapped in PYTHON, and coupled with the affine invariant ensemble sampler (Goodman & Weare 2010) algorithm implemented in the EMCEE package (Foreman-Mackey et al. 2013). The PIKAIA part used 200 individuals (a set of parameters) for 2000 generations, while we run EMCEE with 100 walkers (or chains) for 10000 steps (we remove the initial 2000 steps as burn-in). We repeat this coupled analysis 1000 times.

The same fitting parameters are used in both approaches, that is a linear ephemeris with reference time T_ref and period P_bin, and the LTTE parameters for each k-th body, i.e. a_{k bin}sin i, period P_k, eccentricity e_k, argument of pericentre ω_k, and the time of the passage at pericentre t_{peri, k}. We use the same boundaries of the fitting parameters for this code and the previous one (see Table 5). All the parameters have uniform uninformative priors.

We obtain a large set of solutions, but we select only the solutions that, first, are physically meaningful (i.e. we discard negative eccentricity solutions, since LM is not bounded in the parameter intervals), and have a |$\chi ^2_\mathrm{r} \lt 2$|⁠. For each of these selected simulations, we run a stability⁹ check with rebound and the MEGNO indicator. We run simulations for 10⁵ yr with the whfast integrator and a small step size of 1 d. We apply the full analysis (model fitting with two approaches and stability analysis) and find that all the solutions with |$\chi ^2_\mathrm{r} \lt 2$| are unstable for both data sets.

We show in Fig. 6 the O − C diagram for the two-companion model for the four best solutions (lowest |$\chi _\mathrm{r}^2$|⁠) for both pikaia implementations and both data sets. The four solutions show clearly different contributions from the inner (3) and outer (4) companions, with different periods, amplitudes, and patterns; yet, they fit the observed data points surprisingly well, especially on the ‘reduced’ data set. It is worth noting that both solutions on the full data set are not able to properly reproduce the general trend of the two observing seasons around epoch |$20\, 000$| (1989–1990), being forced to fit the earliest points by Kilkenny et al. (1994).

In Table 6, we present the orbital and physical parameters of these best-fitting solutions. Values for the masses of the companions are within the brown dwarf range. We did not attempt to compute realistic errors (i.e. other than the nominal errors output from the LM fit) on the derived parameters due to the dynamical instability of all the solutions we found.

Additionally, we test a different model with a linear ephemeris (T_c), a one-companion LTTE (τ₃), and a quadratic term (Q). We apply this model to both data sets only with the PIKAIA + EMCEE approach. We use uniform priors within the boundaries in Table 7. We find solutions with |$\chi ^2_\mathrm{r} \gt 6$| (see Table 7 and Fig. 7) and Bayesian Information Criteria that are higher than the two-companion model, for both the data sets. For this reason, we discard this model as a possible explanation for the ETVs.

6 DISCUSSION AND CONCLUSIONS

In this work, we presented a study of the eclipsing binary system HW Vir by using hitherto unpublished photometric observations from four different facilities. We converted all the light curve timings into a common reference frame, as it was crucial for the purposes of this work to have accurate and homogeneous timestamps in order to properly compare different data sets. By combining our new timings with the ones available in the literature, we independently confirmed that the Beuermann et al. (2012) model reproduces the recent literature data until 2011, but it is unable to fit our new timings. Additionally, we tested the dynamical stability of their proposed model and we found it to be unstable after only a few thousand years, opposite to their claim of 10⁸ yr of stability.

As a first effort to find a proper model for the LTTE in HW Vir, we used the PIKAIA code, which implements a genetic algorithm to explore the parameter space and estimate new parameters for the companions of the binary system. We found a set of parameter vectors with a very good fit in a statistical sense, able to explain all the available data. Notwithstanding, these sets of solutions led to very high values for the masses of the companions of HW Vir (∼50M_J, within the mass range of brown dwarfs) and dynamically unstable systems.

Regarding the recent work of Esmer et al. (2021), we describe the most significant differences between their approach and ours in the following. We performed a fully homogeneous analysis of all the new light curves presented, with the same tools and by fitting an accurate EB model (rather than measuring the T₀s with the Kwee & van Woerden 1956 method; Li et al. 2018). This, coupled with the use of larger telescopes, resulted in more accurate eclipse timings by a factor of 5, on average. Also, we exploited a genetic algorithm to perform a comprehensive global search of the parameter space rather than a local one. For this reason, although our search for stable LTTE orbits has been unfruitful, the orbital parameters of our four new solutions fall well outside the region explored by Esmer et al. (2021). The direct O − C comparison of their T₀ with ours is also reassuring, as the average offsets of the residuals measured on a season-by-season basis demonstrate the subsecond accuracy in the absolute timestamp calibration of both data sets.

Although the best-fitting solutions we found were proven to be dynamically unstable, it is worth asking whether other stable orbital solutions with similar LTTE amplitudes exist, and how could we confirm or disprove them with one or more independent techniques.

The prospects for a follow-up with direct imaging are not very promising in the short term. The combination of angular separation (in our best solution, |$0\overset{^{\prime \prime }}{.}11$| and |$0\overset{^{\prime \prime }}{.}47$|⁠, respectively) and contrast (≃10⁻⁵ in the K band if we assume the typical luminosity of a mature 50 M_J brown dwarf; Phillips et al. 2020) falls beyond or very close to the sensitivity limits of the existing ground-based facilities such as SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch; Beuzit et al. 2019) and GPI (Gemini Planet Imager; Ruffio et al. 2017). However, such systems may become very interesting targets for upcoming high-contrast imaging missions such as JWST and the Roman Space Telescope.

On the other hand, astrometry as a follow-up approach could be much more feasible with the release in the near future of the individual astrometric measurements by GAIA (Gaia Collaboration 2016). If we assume that the observed O − C is entirely due to a combination of LTTE signals, its amplitude A_{O − C} can be easily translated into the expected astrometric signal, s, as s = A_{O − C} × c/d, where d is the distance to HW Vir from Table 1. We probe a range of A_{O − C} from 100 to 1500 s, which is spanning the amplitude of the oscillating LTTE terms of the orbital solutions claimed in the recent literature and also compatible with those included in our two best-fitting models in Fig. 6. We find that s ranges from 1.10 ± 0.12 to 16.6 ± 1.8 mas for A_{O − C} = 100 and 1500 s, respectively. That is in principle comfortably within the reach of GAIA sensitivity, since the expected astrometric precision of the individual positional measurements of HW Vir is ∼30 μas (Sahlmann, Triaud & Martin 2015). In such a scenario, the detection will be limited by the temporal baseline rather than the astrometric precision. Yet, if Gaia will survive up to its operational goal of 10 yr, at least the LTTE component with the shortest period can be robustly retrieved, while for the longest one a global analysis combining Gaia with the existing ETV data points will be needed.

A satisfying explanation for the ETVs of HW Vir is still eluding us; however, this only highlights the fact that there is still a lot to be learned about systems of this kind. One of the challenges to accurately determine the underlying cause of the ETVs in this case is that the observations show that the period of one of the components from the LTTE of HW Vir is longer than the total observational time span available. Therefore, increasing the observational baseline will certainly bring us closer to determine the cause behind the ETVs of HW Vir.

ACKNOWLEDGEMENTS

GP and LB acknowledge the funding support from Italian Space Agency (ASI) regulated by ‘Accordo ASI-INAF n. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPS Fasi A/B/C’. LT acknowledges support from Ministero dell'Università e della Ricerca (MIUR; PRIN 2017 grant 20179ZF5KS). DN acknowledges the support from the French Centre National d’Etudes Spatiales (CNES).

DATA AVAILABILITY

The data underlying this article will be uploaded on Vizier/CDS in a second stage; in the meantime, it will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: LITERATURE TIMINGS

In this table, we list the 240 timing measurements taken from the literature (from the compilation by Kilkenny et al. 1994, K94; Lee et al. 2009, L09; Beuermann et al. 2012, B12) and included in our fits together with our new data (Table 3), after being converted by us into a uniform BJD_TDB time standard (Eastman et al. 2010). The epoch is computed according to the ephemeris in equation (1) of (Beuermann et al. 2012).

Author notes