New evidence about HW Vir's circumbinary planets from Hipparcos-Gaia astrometry and a reanalysis of the eclipse timing variations using nested sampling

The post common-envelope eclipsing binary HW Virginis has had many circumbinary companions proposed based on eclipse timing variations. Each proposed solution has lacked in predictability and orbital stability, leaving the origin of the eclipse timing variations an active area of research. Leveraging the catalogue of \textit{Hipparcos} and \textit{Gaia} proper motion anomalies, we show there is slight evidence for a circumbinary companion orbiting HW Vir. We place an upper limit in mass for such a companion which excludes some previously claimed companions. We also apply this method to V471 Tauri and confirm the non-detection of a previously claimed brown dwarf. We adapt the {\tt kima} nested sampling code to analyse eclipse timing variations and re-analyse archival data on HW Vir, varying the order of the ephemeris that we fit for and the amount of the data that we use. Although signals are clearly present, we find two signals around 2500 and 4000 day periods that are not coherent between different \textit{chunks} of the data, so are likely to not be of planetary origin. We analyse the whole dataset and find the best solution to contain four signals. Of these four we argue the outermost is the most compatible with astrometry and thus the most likely to be of planetary nature. We posit the other three pseudo-periodic signals are caused by physical processes on the white dwarf. The eventual release of the full \textit{Gaia} epoch astrometry is a promising way to confirm whether circumbinary planets exist around HW Vir (and other similar systems), and explore white dwarf physics.


INTRODUCTION
Although the majority of known exoplanets have been detected around single stars on the main sequence, planetary systems around post-main sequence stars and in binary star systems are known to exist. The first detected exoplanetary system was around a pulsar (Wolszczan & Frail 1992) and planets orbiting single white dwarfs are known to exist (e.g. Bachelet et al. 2012;Vanderburg et al. 2015Vanderburg et al. , 2020. Many single white dwarf stars have been found to exhibit irregular transit-like and dimming events as well as having atmospheres polluted with heavy elements, both pointing to debris being accreted onto the star which could potentially have been scattered inwards by an invisible companion (Koester et al. 2014;Farihi et al. 2022). Planetary systems around main sequence binaries have been detected in transit by Kepler (e.g. Doyle et al. 2011) and TESS (e.g. Kostov et al. 2020) and also in radial velocity (e.g. Standing et al. 2023).
Planets are therefore known to orbit main-sequence binaries and are able to survive the evolution of a single star. There have been many claims of planets 1 orbiting evolved binaries, but they are yet to ★ E-mail: txb187@student.bham.ac.uk 1 Many of these have masses that would put them above the Deuterium burn-be fully confirmed as planets. These candidate planets, orbiting postcommon envelope binaries, are currently claimed based on periodic variations of the binary's mid-eclipse times. These variations can arise due to the light travel-time effect (LTTE) from the eclipsing binary orbiting the common center-of-mass between itself and the companion. These putative planets could be the counterparts of the detected main-sequence circumbinary planets that have lived through the evolution of their host binary (e.g. Columba et al. 2023).
The existence of these planets has been debated (Mustill et al. 2013). In many cases the claimed planetary solutions fail to predict future eclipses (e.g. Pulley et al. 2022). One candidate, orbiting V471 Tauri (Beavers et al. 1986) was later followed up with direct imaging and was not detected with a high confidence (Hardy et al. 2015).
HW Virginis (HW Vir) is one of the most famous examples of post-common envelope binaries with claimed companions, first proposed by Lee et al. (2009). The system consists of a sdB primary of mass A = 0.418 ± 0.008 M ⊙ and an M-dwarf secondary of mass B = 0.128 ± 0.004 M ⊙ in a binary with orbital period bin = 0.116719556 ± 7.4 × 10 −9 days 2 . Eclipses have been precisely measured for over 30 years, with many conflicting solutions proposed (e.g. Beuermann et al. 2012;Esmer et al. 2021) with either one or two planets proposed as the cause of the eclipsing timing variations. One major issue is that none of the single-planet solutions fit the data satisfactorily, but none of the better-fitting two-planet solutions appear to be dynamically stable (Brown-Sevilla et al. 2021;Mai & Mutel 2022). Another issue, as mentioned above, is that all of the proposed solutions very quickly diverge from the data subsequently collected.
Non-planetary explanations have been suggested which can produce eclipse timing variations in short-period binaries such as HW Vir. The period (or apparent period) of the binary could be affected by apsidal precession if it is eccentric, and magnetic braking (Rappaport et al. 1983) or emission of gravitational waves (Paczyński 1967) could cause the orbit to shrink due to angular momentum loss. Other magnetic effects have also been proposed, such as the Applegate mechanism (Applegate 1992), or a more recent mechanism, requiring less energy suggested by Lanza (2020). However in most cases these are insufficient to fully explain the shape or the amplitude of the observed modulations in eclipse time.
Many of these candidate planets will need to be confirmed/rejected through other methods. One example of this happening is with V471 Tau. This system consists of a WD primary of mass A = 0.797 ± 0.016 M ⊙ and a K-Dwarf secondary of mass B = 0.864±0.029 M ⊙ in a binary with orbital period bin = 0.5211834194±7.2×10 −9 days (Muirhead et al. 2022). This system shows periodic variations of the mid-eclipse times, which have been used to suggest an orbiting brown dwarf (Beavers et al. 1986;Guinan & Ribas 2001). The system has since been directly imaged with SPHERE, and these observations resulted in a non-detection (Hardy et al. 2015), thus rejecting the claimed brown dwarf.
Planets around ultra-short period evolved binaries such as these may also eventually be detectable in gravitational waves, for example by the Laser Interferometer Space Antenna, LISA (Danielski et al. 2019). LISA will however only be sensitive to binaries of shorter orbital period than HW Vir.
Another possibility for confirming or rejecting post-common envelope circumbinary planets is precise astrometry. The space telescope Gaia (Gaia Collaboration et al. 2021) is performing a precise astrometric survey of the whole sky which will have a baseline of about 10 years. Gaia's astrometric solution will be able to investigate some of these systems without relying on any eclipse timing data (Sahlmann et al. 2015). However, individual astrometric measurements are expected to released around late 2025. In the meantime, we can only rely on the proper motion anomaly method (Kervella et al. 2019). Before Gaia, Hipparcos (ESA 1997) also performed an astrometric survey, but of a much smaller sample of stars, at a lower precision. HW Vir is within that sample, as is V471 Tau. The proper motion anomaly method combines positions and proper motions of a star from Hipparcos and one of Gaia's recent data releases to estimate the effect of an orbiting of an orbiting companion on the proper motion of the star. This method has been applied to single stars, and combined with other techniques such as radial velocity, has led to the detection and characterisation of several planetary companions (e.g. Mesa et al. 2022;Rickman et al. 2022).
In this paper we present a new piece of astrometric information, in the form of the proper motion anomaly, to the puzzle that is HW Vir, 2 Parameters are taken from Esmer et al. (2021), these values are used for the rest of the analysis when the mass of the central binary is needed and we perform a new fit of the eclipse timing data utilising nested sampling and analysing different chunks of data separately. We report on the lack of consistency of signals in the eclipse times and present the model that we find to fit best to the whole dataset, providing a suggestion of which signal is favoured to be planetary since not all the detected signals can be.
The paper is set out as follows. We describe the proper motion anomaly method, and apply it to both V471 Tau and HW Vir in Section 2. In Section 3, the use of kima to fit eclipse times is described. Section 4 details the eclipse timing data used, and the results from the analysis of the data. We discuss the results and implications and conclude in Section 5.

ASTROMETRY: USING THE PROPER MOTION ANOMALY METHOD
We firstly explain the method of the astrometric proper motion anomaly, and secondly apply this to both V471 Tau and HW Vir. We compare the results from the astrometric proper motion anomaly to some of the previously proposed planetary solutions.

How does the proper motion anomaly work?
The proper motion anomaly analysis method is described in detail in Kervella et al. (2019). Using positional measurements from Gaia and Hipparcos, we determine the long-term, mean proper motion vector of the system HG by dividing the observed change in position by the time baseline HG between the two measurements (that is, 24.75 years between Hipparcos and Gaia DR3). For the nearest stars, second order effects must be taken into account in this computation, but they are negligible for the systems discussed in the present paper. Thanks to the long time baseline, and assuming that the orbital period of the companion is significantly shorter than HG , HG essentially traces the proper motion of the center of mass of the system. Separately, the short-term proper motion measurements Hip and DR3 (obtained respectively by Hipparcos and Gaia) trace the vector sum of 1) the linear proper motion of the center of mass and 2) the orbital motion orbit of the photocentre of the system around the center of mass. Figure 1 visually presents the different vector quantities considered in these computations.
When considering a planetary companion, the photocentre is located very close to the geometrical center of the star. In this configuration, subtracting the long-term proper motion HG from the Gaia DR3 short-term proper motion DR3 gives access to the proper motion anomaly of the star Δ , that traces the orbital motion of the star around the center of mass of the system. The quantity Δ can be scaled to a linear tangential velocity anomaly tan using the parallax. This is the two-dimensional counterpart of the radial velocity that is traditionally employed to detect exoplanets.
Corrective terms must be considered to interpret the measured proper motion anomaly in terms of companion properties. Firstly, the orbit quantity is an average quantity over the Gaia integration window, that has a duration of DR3 = 34 months. This averaging implies that the measured proper motion anomaly will be smeared, reducing the sensitivity in terms of companion mass. This loss of sensitivity is particularly strong for orbital periods shorter than DR3 Secondly, the time baseline HG between Hipparcos and Gaia DR3 results in the subtraction of part of the proper motion signature of very long period companions ( > 3 × HG ) during the computation of the HG quantity. This effect induces a loss of sensitivity to such very long period companions. These two effects determine the companion mass sensitivity function of the proper motion anomaly method.
We now derive the equation for the tangential velocity caused by a companion if measured instantaneously. This equation can then be combined with the sensitivity function calculated numerically. The proper motion is usually divided into its components in rightascension (ra) and declination (dec).
= ra e ra + dec e dec , with e ra and e dec the basis vectors in ra and dec. We then subtract the long-term HG proper motion from the Gaia DR3 proper motion and take the magnitude of this vector to get the tangential velocity anomaly.
= ( DR3,ra − HG,ra )e ra + ( DR3,dec − HG,dec )e dec , where is the parallax. Now given an inner mass and an outer mass the relative orbital velocity is where G is the gravitational constant, the semi-major axis of the relative orbit, and the relative orbital distance at the measured time. The distance is given by with and the eccentricity and true anomaly of the orbit at the measured time (they are the same for the relative orbit or the orbit of one of the components). Combining these two equations and using that the velocity of the inner body (i.e. the luminous one) relates to the outer velocity by gives us: where 1 is the semi-major axis of the outer orbit which relates to that of the relative orbit by 1 = + .
This derivation is valid for an instantaneous measurement of tan , which would correspond to an instantaneous measurement of DR3 and an infinitely long baseline for HG . This is, of course, not the case. As described above we must also include a sensitivity function. This function has been numerically calculated by Kervella et al. (2019). This leads to the sensitivity curve for the proper motion anomaly at different periods which is used in the following section. These curves give the areas of Period-Mass space that are consistent with a measured proper motion anomaly, under the assumption of a circular orbit.

Applying the proper motion anomaly to V471 Tau and HW Vir
Calculating the long-term proper motion between Hipparcos and Gaia DR3 has been done and combined into a catalogue by both Kervella et al. (2022) and Brandt (2021), using different combinations of the main Hipparcos reductions. The values for the tangential velocity for HW Vir and V471 Tau are shown in Table 1. We note that the values are in good agreement between both catalogues and choose arbitrarily to use the Kervella et al. (2022) value from now on.

V471 Tau
For V471 Tau, we have a proper motion anomaly between 2-and 3-. Figure 2 shows the sensitivity curve of the proper motion anomaly method associated with the value for V471 Tau from Kervella et al. (2022). The dark green line shows the curve on which a body needs to lie to produce the observed tangential velocity value, the darker and lighter shaded regions show the 1 and 3 regions around that line. The spikes towards shorter orbital periods are the result of the sensitivity function as described above. The slope at longer orbital periods is produced when only small fraction of an orbital arc is covered and hence the efficiency function is small. In between, is a region of highest sensitivity. This provides an upper bound that is far below the mass of the proposed solutions by Beavers et al. (1986) and Guinan & Ribas (2001). This re-affirms the conclusion of Hardy et al. (2015) which did not find evidence of the proposed brown dwarf, and confirms that the variations in the mid-eclipse times must be coming from some other source.
We numerically estimate the proper motion anomaly that would be caused by the binary. For a given set of parameters ( 0 , 1 , ) we perform a bisection algorithm suggesting values for tan , and comparing the value of 1 obtained (given 0 and ) to the given value, until the masses agree to 0.001M jup . We repeat this for 1000 realisations of the binary parameters to then obtain the median and 1 values of 32 +13 −21 m s −1 . This is entirely consistent with the tentative signal that is seen. The proper motion anomaly is sensing the smeared binary motion and so does not suggest an orbiting companion.

HW Vir
For HW Vir, the tangential velocity is distinct from zero at around 2 confidence in both catalogues. We cannot therefore conclude from the proper motion anomaly that there is definitely an orbiting body, but this Gaia-Hipparcos combined measurement brings new evidence that suggests such a body is more likely to exist than not.
First, we validate that this tentative proper motion anomaly is not caused by the smeared orbital motion of binary. In the same way as for V471 Tau, we numerically estimate the proper motion anomaly  that would be induced by the binary and obtain the median and 1 values of 2.52 +0.80 −1.60 m s −1 . The excess tangential velocity is therefore not caused by the HW Vir binary.
The top panel of Figure 3 shows the sensitivity curve of the proper motion anomaly method associated with HW Vir. The curve has the same shape as in Figure 2 (since all the spikes are primarily related with the Gaia 34-month observing window) but is zoomed in on the area of best sensitivity. We overplot the locations of the orbiting bodies proposed by three previous studies (Beuermann et al. 2012;Brown-Sevilla et al. 2021;Esmer et al. 2021). We note that four of the proposed solutions include one orbiting body above the 3 line. These solutions are disfavoured 3 by the observed proper motion anomaly which is too weak to have been produced by these putative objects. This plot also shows the locations of the four components from our best-fitting model 4 , which we describe in Section 4.2.2.
This tentative proper motion anomaly is an extra piece of information about the HW Vir system which provides astrometric evidence that there may be an orbiting circumbinary companion.
The catalogues of accelerations mentioned above rely on Gaia EDR3 (Gaia Collaboration et al. 2021). The same analysis was done earlier using Gaia DR2 (Gaia Collaboration et al. 2018), and the value from Kervella et al. (2019) for HW Vir is 309 ± 200. From this 3 They may still be possible in reality, if we have a very eccentric orbit for the companion, and we observe it close to apastron (where the motion is slower) 4 We do not claim that all four of the signals are indeed planets we infer that the astrometric signal is getting more confident as more Gaia data becomes available. This implies that if there is indeed a signal there, it should be detectable from future Gaia data releases.

FITTING ECLIPSE TIMING VARIATIONS WITH KIMA
Whilst verifying whether proposed solutions for the HW Vir systems were compatible with the proper motion anomaly, we also decided to re-analyse the eclipse times of HW Vir with a nested sampler, which we believe has not been attempted yet. Most of the literature uses 2 maps or reduced 2 to make inferences about the number of signals present in the data, but none have conducted a Bayesian model comparison in this way yet.
Amongst Bayesian methods, nested sampling has the advantage to let some key parameters free that are usually fixed in other types of analyses. In our case, the number of orbiting planets, p is a free parameter which allows for a robust model comparison, based on a ratio of Bayesian evidence. All planetary signals are adjusted at once and models with 0, 1, 2... planets are constantly compared to one another.
kima is an orbital fitting algorithm originally designed for application to radial velocities (Faria et al. 2018). We adapt it to fit mid-eclipse times instead, to then apply it to HW Vir. kima leverages nested sampling using DNEST4 (Brewer & Foreman-Mackey 2018) to explore parameter space and calculate the likelihood of proposed samples. Using the trans-dimensional sampling in kima the number of Keplerian signals, p , being fit is a free parameter as described above. This allows a comparison of the different numbers of signals present in the data with a Bayes Factor. The Bayes Factor for a p = model 5 compared to a p = − 1 model is the ratio of the evidence for each model. The evidence is the primary output of nested sampling and is the integral of the likelihood over the prior mass. In nested sampling this integral is calculated as a weighted sum, with the weights being associated to the change in prior mass between consecutive samples (Skilling 2006). In this case the evidence for an p = model is the sum of the weights of all the samples with planets, and then the Bayes Factor is = −1 . We use a detection threshold of 150 as recommended by Kass & Raftery (1995). A Bayes Factor larger than 150 is taken as very strong evidence for the more complex model over the less complex model (and roughly equivalent to a p-value of 0.001). It is common that a nested sampler finds the sum of the weights of all the samples is highest for the highest p explored by the sampler (Faria et al. 2018;Standing et al. 2022). However, so long as the ratio is not > 150 those most complex solutions, whilst providing a better fit to the data, do not contain enough statistical evidence to warrant the extra number of parameters.
As a by-product of the nested sampling to calculate the evidences, we can obtain posterior samples for the various parameters from kima. These allow us to perform parameter estimation on any detected signals, assuming a light travel-time effect (LTTE) model with a companion on a Keplerian orbit. kima has already been used to detect and test the detectability of circumbinary planets with radialvelocities (Triaud et al. 2022;Standing et al. 2022Standing et al. , 2023. We redirect the reader to these publications for more thorough explanations of how the model comparison works.
We fit the eclipse times in kima with a number of Keplerian signals 5 p meaning number of planets in the model as well as an ephemeris function. We allow the ability to fit for one of a linear, quadratic or cubic ephemeris. These are shown in equation 9 below: where is the epoch of an eclipse (i.e. the number of the eclipse with the first eclipse being 0), ( ) is the time of that eclipse in our model, 0 is the reference time (time at epoch 0 here), 0 is the period of the eclipsing binary at the reference time, 0 and 0 are the first and second time-derivatives of the binary period (at the reference time), and is the time-delay due to the LTTE of an orbiting body. The middle three terms are a Taylor series and if we ignore the terms of order ≥ 2 we are using a linear ephemeris, if we ignore terms of order ≥ 3 we are using a quadratic ephemeris, and using all the terms above is a cubic ephemeris. The functional form for the LTTE due to an orbiting body is as in Irwin (1952): where and are the eccentricity and argument of periastron of the orbiting body, ( ) its true anomaly at time and the semiamplitude of the signal: where and are the mass and orbital period of the orbiting body, the total mass of the eclipsing binary, i the inclination to the line of sight of the planetary orbit, and and the speed of light and gravitational constant.
In equation 10, and are functions of (time). The orbital period of an orbiting body is much greater than the difference in time from all other terms, so we use ≈ 0 as a first order approximation.

FITTING ECLIPSE TIMES
In this section we first detail from where the eclipse timing data is obtained, and then present the results from the analysis using kima.

Data for HW Vir
We use archival data for HW Vir eclipse times (considering only the primary eclipse). We use the data from Brown-Sevilla et al. (2021), which collated data from Kilkenny et al. (1994), Lee et al. (2009), andBeuermann et al. (2012), as well as their own data. To this we add the data from Baran et al. (2018), Esmer et al. (2021), and Mai & Mutel (2022). Of the data reported in Baran et al. (2018), we find that the data taken with SAAO have a small offset of ∼ 80 sec from the rest of the datasets (including the other data reported in the same publication). Since there is still good coverage without this, we exclude these data from our analysis.
We perform the analysis on the whole dataset, but also divide it into smaller chunks to assess how consistent any signals that appear are. This way we can assess if, although the overall model doesn't have good predictive power, a subset of the signals might be predictably and consistently present. We divide the dataset into chunks of approximately 1/3 and 2/3 the length of the whole dataset, with epochs as shown in Table 2. The chunk "tier3", is extended back an extra 10 000 epochs and overlaps with "tier2". The different chunks are also visualised alongside the data in the top panel of Figure 4.

Results from eclipse timing variation fits
In this section we present the results from a reanalysis of the mideclipse time data, analysed using kima. The nested sampling implemented requires a prior distribution for each parameter, these are detailed in Table 3. The Kumaraswamy distribution (Kumaraswamy 1980) approximates the beta distribution, and the shape parameters as shown in Table 3 are those that Kipping (2013) argues best represent the distribution of exoplanetary eccentricities, based on exoplanets detected with the radial velocity method. The analysis is performed with each of a linear, quadratic, and cubic ephemeris for each chunk as well as for the whole dataset.

Results from analysing different chunks
A LTTE signal due to an orbiting companion should to be coherent in time. The analysis of different chunks, using a Keplerian prescription, would therefore be expected to lead to posteriors that are consistent across chunks. Our analysis of the different chunks shows a lack of consistency and therefore casts doubts on the ETV signals being solely due to an orbiting companion (or more). Throughout the analysis of the different chunks of data, recurring signals are seen around two periods: 4000 days and 2500 days. Longer period (and higher amplitude) signals do exist in many of the chunks, however they are far from consistent.
To assess the consistency of the signals at the recurring periods, the clustering algorithm HDBSCAN (McInnes et al. 2017) is used to identify clusters in the P-K plane from the kima posterior samples. The clusters are then visually associated with one of the two recurring periods or not. The lower panels of Figure 4 show the clusters of posterior density around 2500 and around 4000 days from runs of kima on different chunks of data 6 . These are shown as corner plots (Foreman-Mackey 2016) between the period and semi-amplitude . While there is a cluster of posterior density around this period in each 7 of these runs, the periods and amplitudes of the signals vary between the runs.
The lack of consistency of these signals points to them not being due to a Keplerian LTTE orbit. These signals may then have a non-6 tier3 did not show signature of a detectable signal around 2500 days so only 5 chunks are shown 7 tier3 notwithstanding periodic or quasi-periodic source. If this is the case, then attempting to fit them with strictly periodic Keplerian signals is unideal. This is exemplified in the upper panel of Figure 4, where we show the best fitting model from the run where tier1-3 is analysed using a linear ephemeris. The best model, a sum of three Keplerian orbits, is woefully incorrect for the middle section. This shows how not only are Keplerian LTTE models not successful at predicting future eclipse times, but they are unsuccessful at interpolating.
In the future, a better approach might be to use Gaussian Processes (Rasmussen & Williams 2006) to model the shorter eclipse timing variation signals. These tools are particularly good at modelling quasi-periodic functions to stellar activity for instance (in photometry and spectroscopy; Barros et al. 2020;Faria et al. 2016) and would seem appropriate in the case of HW Vir.

Results from analysing the full dataset
We now show the results from an analysis of the full dataset. We allow kima to fit freely up to p = 6 signals along with a quadratic ephemeris. One advantage of using kima is its ability to assess the number of signals present using Bayesian model comparison. In this case a four-signal solution is favoured as it is the highest number of planets with a significant Bayes Factor over a model with one fewer planet. The respective Bayes Factors can be seen in Table 4. Four signals is more than most other analyses which only find up to two signals. The two signals already discussed (around 2500 and 4000 days) are both present in the best-fitting solution. We know this because a large fraction of the posterior sample congregate at these two orbital periods (as shown in the upper panel of Figure 3). The other two signals are not nearly as well constrained and do not correspond to any clear over-density in the posterior, likely because these are longer signals that have not had the chance to repeat yet, making their parameters uncertain.
Past analyses have regularly identified a signal corresponding to the one we find around 4000 days (e.g. Beuermann et al. 2012;Esmer et al. 2021), however none have identified a signal near 2500 days.
We do not consider any formal stability arguments. Many previous studies have found that multiple-planet solutions are unstable, and since we have a strong reason to doubt that the detected signals Top-right: all four signals included, the most massive is shown in green and the sum of all four in blue, this signal in green is the one we claim as the most likely candidate for being a planet. The x-axis is the Epoch as in equation 9. Bottom-right: the inner three Keplerian functions are shown with the sum of these three in blue. The data is represented with the fourth, large-amplitude, signal removed. These three signals are most likely not of planetary nature. within the eclipsing times are produced by an orbiting planet, we feel a stability analysis is meaningless. The upper-left diagram of Figure  5 shows the orbital configuration of planets corresponding to all four signals. The inner two of them are reasonably circular, the outer two are more eccentric, and the outermost crosses the other orbits.
Clearly not all four signals can be from orbiting bodies. Ignoring the outermost, eccentric orbit, the lower-left diagram in Figure 5 shows the orbital configuration of planets corresponding to the inner three signals. While these three signals do not cross into each others orbits, they present a very compact configuration, that would likely not be stable either.
The astrometric tangential velocity implies it is more likely than not there is one orbiting companion to the HW Vir binary. Of the four signals, if one is of planetary nature, we favour the fourth and outermost signal. The analysis reported in section 4.2.1 casts strong doubts on the inner two signals since they appear only quasi-periodic. The third signal has too long a period to assess the consistency with the chunking method, but its 'orbital parameters' are similar to the inner two signals with a small amplitude and mild eccentricity. Compared to all others, the outermost signal lies closest to where the median value of the proper motion anomaly predicts (the dark line on Figure 3). We note that this candidate planet signal is of a similar mass and period to components of the solutions by Esmer et al. (2021) and Brown-Sevilla et al. (2021). These all likely correspond to the same signal but vary in orbital period due to the data having not covered multiple cycles yet.
The plots on the right hand side of Figure 5 show the model curves for each of the signals along with the combined model and data. The Table 5. Parameters from the analysis of the full dataset, with a quadratic ephemeris. The keplerian parameters for each signal is shown as if it was keplerian LTTE orbit. a for the circular orbit we combine the and parameters together since otherwise they are extremely correlated and no information can be gained. b for the two outer signals, the posterior density is not well constrained so clusters around the best-fitting solutions are used for some of the parameters. We note the uncertainties are likely too small to represent the true uncertainty in the model. c For these periods and amplitudes, the uncertainty is reported as the difference between the value of the parameter when a linear ephemeris is fit and when a quadratic ephemeris is fit. The median of the posterior density cluster from the quadratic ephemeris fit is retained as the quoted value. The mass distribution is the propagated in a monte-carlo way.

Parameter
Value Units upper-right panel shows the full model in the background and the outer orbit Keplerian signal in green, the lower-right panel shows the other three individual signals in shades of purple as well as their sum in the background. The three signals shown together in the lower right panel are those we claim to be most likely not produced by a planet (especially the two at shorter periods), these might be better modelled together as a Gaussian Process.
While we know that this four-component solution cannot correspond to four orbiting companions, to allow future comparison with our work, we still report the parameters of the orbits as if they were real. The parameters are detailed in Table 5. The uncertainties asso-ciated with the parameters of the inner two orbits are well-defined as they are associated with clear clusters of posterior density (clustering using HDBSCAN is also used here). The outer two orbits do not belong to large clusters, so while they can be associated with clusters found by HDBSCAN the uncertainty on the parameters extracted from these are likely underestimated. This is because there has not been enough data for the signal to repeat. In the case of the outer signal the data has not even covered a whole phase yet. This also causes a degeneracy between the orbital parameters of the outer signals and the ephemeris terms. To partly address the underestimation, we take two analysis runs with kima, one using a linear ephemeris and one using a quadratic. The uncertainty on the amplitude, period are then taken as the difference between the values from the two models, with the quoted value remaining the value from the analysis with a quadratic ephemeris. This is to keep the whole table representing a coherent solution. Corresponding planetary masses for the outer two signals are then produced in a monte-carlo way. It should be noted this is therefore not a statistically derived uncertainty, but is a rough representation of the uncertainty from the fitting procedure.

DISCUSSION AND CONCLUSION
Our analysis of the eclipse times of HW Vir does not find a single conclusive solution. This is in agreement with past work since every published solution has subsequently diverged from new data acquired afterwards (e.g. Pulley et al. 2022). We have shown that there are two strong periodicities in the full dataset which are also seen independently in some of the smaller chunks. However, though signals can be found near these periods in most of the chunks, their posterior distributions in and are not completely coherent in time, nor statistically consistent with one another. While past analyses have identified the periodicity around 4000 days, none have identified the signal around 2500 days.
We have presented our solution from the analysis of the whole eclipse timing dataset performed in a fully Bayesian way using nested sampling within kima. This 4-component model includes signals at both of the strong periodicities. It is abundantly clear that not all four of the signals in this model are due to orbiting companions, in fact it is possible that none of them are. There likely must be some other mechanism involved for the variation in eclipse times, one possibility being a magnetic effect which is not yet fully understood (e.g. Lanza 2020). We suggest that using a Keplerian prescription for non-planetary, quasi-periodic signals like what these appear to be is insufficient and that using a Gaussian Process method may work better (as in Faria et al. 2016). We propose that of the four signals, if one is produced by a planet, it is most likely to be the outermost one. This signal also best fits the astrometric evidence and has a signature that looks most different to the other three.
We have applied the proper motion anomaly method to V471 Tau and confirmed the non-detection of a previously proposed orbiting brown dwarf. We have also applied it to HW Vir and shown that there is a tentative 2-signal of an acceleration due to an orbiting body. From the upper limit this poses, we can discount some previously proposed companions which are too massive to be consistent with the proper motion anomaly. Comparing the astrometric signal with the four signals we extract from the eclipse timings in Fig. 3, we find that the outermost signal is the most consistent. If correct this corresponds to a 17 M jup , 16 000 day, highly eccentric companion. Thanks to additional data, a longer baseline, and an improved astrometric solution, the full epoch astrometry from Gaia will (circa 2025) likely be able to help resolve whether the HW Vir binary is indeed host to an orbiting circumbinary companion.The Gaia baseline will still be much shorter than the most likely planet's orbital period, so while the whole period would not be covered, astrometry may still tell us whether or not such a planet exists, independently of the ETVs. This will help identify which (if any) of the varying eclipse timing signals is actually caused by that orbiting body. In turn, this will help isolate the functional form of the potential new physics causing the other signals (for instance the 2500 and 4000 day signals). As described in Sahlmann et al. (2015), thanks to Gaia's final solution, other post-common envelope circumbinary systems will be solved astrometrically with our paper being the first attempt at doing so.