Transit timing variations in the WASP-4 planetary system

TL;DR

This paper confirms transit timing variations in WASP-4, suggesting possible orbital decay, and rules out common causes, contributing to understanding hot Jupiter dynamics.

Contribution

It provides the first high-confidence detection of TTVs in WASP-4 and discusses potential causes like orbital decay and precession.

Findings

Transit times are 80s earlier than expected.

Confirmed TTVs fit a quadratic ephemeris with period decay.

Orbital decay and precession are plausible causes.

Abstract

Transits in the planetary system WASP-4 were recently found to occur 80s earlier than expected in observations from the TESS satellite. We present 22 new times of mid-transit that confirm the existence of transit timing variations, and are well fitted by a quadratic ephemeris with period decay dP/dt = -9.2 +/- 1.1 ms/yr. We rule out instrumental issues, stellar activity and the Applegate mechanism as possible causes. The light-time effect is also not favoured due to the non-detection of changes in the systemic velocity. Orbital decay and apsidal precession are plausible but unproven. WASP-4b is only the third hot Jupiter known to show transit timing variations to high confidence. We discuss a variety of observations of this and other planetary systems that would be useful in improving our understanding of WASP-4 in particular and orbital decay in general.

Tables

Table 1. Table 1: Times of mid-transit for WASP-4 determined in this work. The final column shows ‘y’ if spot crossings are clearly visible in the data.

Telescope	Date	Epoch	BJD(TDB)		Spot
Danish	2008/08/19	$-$1351	2454697.79821 $\pm$	0.00010
Danish	2008/08/22	$-$1348	2454701.81293 $\pm$	0.00011	y
Danish	2008/09/22	$-$1325	2454732.59192 $\pm$	0.00017
Danish	2008/10/01	$-$1319	2454740.62150 $\pm$	0.00007	y
Danish	2009/08/25	$-$1073	2455069.82652 $\pm$	0.00009
Danish	2011/08/03	$-$544	2455777.75045 $\pm$	0.00014	y
Danish	2011/08/11	$-$538	2455785.78063 $\pm$	0.00011	y
NTT	2011/10/21	$-$485	2455856.70662 $\pm$	0.00007
Danish	2015/08/22	562	2457257.83498 $\pm$	0.00007
Danish	2015/08/26	565	2457261.84940 $\pm$	0.00010
Danish	2015/08/30	568	2457265.86449 $\pm$	0.00010
Danish	2016/08/12	828	2457613.80465 $\pm$	0.00010
Danish	2016/08/20	834	2457621.83379 $\pm$	0.00011
Danish	2016/08/24	837	2457625.84865 $\pm$	0.00010	y
SAAO 1.0m	2016/10/13	874	2457675.36302 $\pm$	0.00016
SAAO 1.0m	2016/10/17	877	2457679.37808 $\pm$	0.00018
Danish	2017/08/07	1097	2457973.78903 $\pm$	0.00018
Danish	2017/08/27	1112	2457993.86231 $\pm$	0.00014
Danish	2018/08/14	1375	2458345.81705 $\pm$	0.00012
Danish	2018/08/18	1378	2458349.83148 $\pm$	0.00010	y
Danish	2018/08/26	1384	2458357.86123 $\pm$	0.00008
SAAO 1.0m	2019/06/18	1605	2458653.61010 $\pm$	0.00020

Equations

$$T_0=\mathrm{BJD}(\mathrm{TDB})2456505.748953(20)+1.338231429(20)E$$

$$T_0$$

$$\frac{dP}{dt}=-9.2\pm 1.1\text{mbox}{msyr}^{-1}$$

Full text

1 Introduction

Thousands of transiting planets are currently known. Under the assumption of Keplerian motion, the transits are expected to recur with a strict periodicity. However, there is a range of phenomena that might change this behaviour, including stellar activity (Oshagh et al., 2013), gravitational interactions (Holman & Murray, 2005), the light-time effect in multiple systems (Woltjer, 1922), and orbital decay due to tides (Birkby et al., 2014). Because of this, extensive efforts have been devoted to the detection of transit timing variations (TTVs) (e.g. Gibson et al., 2010; Harpsøe et al., 2013; Maciejewski et al., 2018b). Many examples have been found for small planets observed using the Kepler satellite (Mazeh et al., 2013) but only two detections have been made for a planet outside the Kepler field. Maciejewski et al. (2016) used transit timing measurements taken over 3.3 yr to find a quadratic term in the orbital ephemeris of WASP-12 to a significance of 5$\sigma$ (see also Patra et al., 2017; Maciejewski et al., 2018b), and WASP-47 b shows sinusoidal TTVs due to the presence of two smaller short-period planets (Becker et al., 2015).

In this work we show that WASP-4 b (Wilson et al., 2008) is the third giant planet known to exhibit detectable TTVs, with a significance of 8.4$\sigma$. This confirms preliminary suggestions from Bouma et al. (2019, hereafter B19) based on data from the TESS satellite (Ricker et al., 2015). Baluev et al. (2015) have also claimed a tentative detection of TTVs in WASP-4, in the form of a sinusoidal variation with a period near 5 d; we do not confirm this detection. When needed, we will adopt the physical properties of the WASP-4 system given by Southworth (2012). In brief, the planet has mass 1.25$\,{\rm M}_{\rm Jup}$ and radius 1.36$\,{\rm R}_{\rm Jup}$, and is on a 1.338 d orbit around a star of mass 0.93$\,{\rm M}_{\odot}$ and radius 0.91$\,{\rm R}_{\odot}$. The star’s effective temperature is 5500 K (Mortier et al., 2013) and it shows magnetic activity in the form of dark starspots (Sanchis-Ojeda et al., 2011). No companion stars have been found in high-resolution imaging studies (Ngo et al., 2015; Evans et al., 2016) and no long-term trend in systemic velocity has been identified (Knutson et al., 2014).

2 Observations

We have been monitoring transits in the WASP-4 system for approximately a decade using the Danish 1.54 m telescope and DFOSC imager at ESO La Silla, as a side-project of the MiNDSTEp Consortium (Dominik et al., 2010). Four transits observed in the 2008 season were presented in Southworth et al. (2009b), and we also observed one transit in 2009, two in 2011, three in 2015, three in 2016, two in 2017 and three in 2018. All observations were obtained through a Cousins $R$ filter with the CCD unbinned and the telescope significantly defocussed to enhance the photometric precision (Southworth et al., 2009a). The scatters of the datasets versus transit fits (see below) range from 0.54 to 1.03 mmag with cadences of typically 140 s.

We also observed one transit in the WASP-4 system using the NTT and EFOSC2 imager (Buzzoni et al., 1984) with an SDSS $r$ filter (ESO filter #784) on the night of 2011/10/21. The telescope was defocussed, the CCD was not binned, and only one comparison star was used due to the small field of view. The light curve is of very high quality (see Tregloan-Reed & Southworth, 2013), with a scatter of 0.47 mmag.

Finally, we observed three transits using the SAAO 1.0 m telescope at Sutherland, South Africa: two in 2016 and one in 2019. The STE4 CCD was used with a Cousins $R$ filter, 2$\times$2 binning to decrease the readout time, and the telescope defocussed. The resulting light curves have scatters of 1.5, 2.3 and 1.6 mmag. In the most recent run (2019 June) we cross-checked timings with another telescope at the observatory to confirm their reliability.

2.1 Data reduction

The raw data were reduced into light curves using the defot pipeline (Southworth et al., 2009a, 2014) which implements aperture photometry, image motion tracking by cross-correlation with a reference image, and calculation of a differential-magnitude light curve. The light curve was obtained by simultaneously fitting the coefficients of a low-order polynomial versus time, and the weights of a set of comparison stars, to the data obtained outside transit. In all cases we used a linear or quadratic polynomial and the coefficients were included as fitted parameters in all subsequent analyses.

The best light curve was obtained by manually iterating the aperture sizes and which comparison stars were included – we found that the choices made here do affect the scatter of the light curve but have no significant effect on its shape. We did not perform bias or flat-field calibrations because they have a negligible effect on DFOSC data except for an increase in the noise level of the light curve (Southworth et al., 2014). All new data are shown in Fig. 1.

The timestamps were converted to the BJD(TDB) timescale using routines from Eastman et al. (2010). Manual time checks were performed for images in most observing sequences, and no discrepancies were noticed111The 2009 season with the Danish Telescope suffered from a timing problem (Southworth et al., 2009c, 2010). We have investigated this further and found that it did not affect any observations of WASP-4 presented here.. We have previously found a good agreement between the timestamps on the Danish and MPG 2.2m telescopes (Southworth et al., 2015), supporting the reliability of both.

Southworth et al. (2009b) used an early version of defot so the 2008 data were re-reduced with the current version. This allowed the use of image motion tracking, additional comparison stars, and timestamps on the BJD(TDB) timescale. The re-reduced data supersede those from Southworth et al. (2009b).

3 Measurements of the transit times

Each transit was fitted in isolation using the jktebop222jktebop is written in fortran77 and the source code is available at http://www.astro.keele.ac.uk/jkt/codes/jktebop.html code (Southworth, 2013, and references therein). The fitted parameters were the transit midpoint $T_{0}$, the fractional radii of the two components ($r_{\rm A}=\frac{R_{\rm A}}{a}$ and $r_{\rm b}=\frac{R_{\rm b}}{a}$ where $R_{\rm A}$ is the radius of the star, $R_{\rm b}$ is the radius of the planet, and $a$ is the semimajor axis of the relative orbit) expressed as their sum and ratio, the orbital inclination, and the coefficients of the polynomial versus time. We modelled limb darkening using the quadratic law, fitted for the linear coefficient, and fixed the quadratic coefficient at an appropriate value.

Uncertainties on the transit times were obtained using Monte Carlo simulations and multiplied by $\sqrt{\chi_{\nu}^{\,2}}$ where $\chi^{2}_{\nu}$ is the reduced $\chi^{2}$ of the fit. This final step is necessary because the aperture photometry routine used in our pipeline tends to underestimate the observational uncertainties. Table 1 gives the 22 transit times calculated in this work. To this dataset we added all timings from table 2 of B19, which includes 41 measurements from the literature and 18 timings from the TESS light curve.

4 Transit timing analysis

Our analysis was performed on the transit timings assembled above. We fitted a straight line to the times as a function of orbital cycle, using a zeropoint for cycle number near to the midpoint of the available data to avoid correlations between the parameters of the ephemeris (Fig. 2). This yielded a poor fit so we tried a quadratic ephemeris instead, finding a significant improvement. The best-fitting linear ephemeris is:

[TABLE]

with $\chi_{\nu}^{\,2}=1.80$ and an rms scatter of 24.6 s. The bracketed quantities indicate the uncertainties in the preceding digits and have been increased by $\sqrt{\chi_{\nu}^{\,2}}$ to account for the imperfect fit. The best-fitting quadratic ephemeris is:

[TABLE]

with $\chi_{\nu}^{\,2}=1.41$ and an rms scatter of 19.9 s. As above, the uncertainties have been increased by $\sqrt{\chi_{\nu}^{\,2}}$ to account for the imperfect fit. We also fitted a cubic ephemeris but found it to be a negligible improvement on the quadratic ephemeris.

This quadratic coefficient corresponds to a period change of $\frac{dP}{dE}~{}=~{}(-3.91\pm 0.46)\times 10^{-10}$ day per orbital cycle, and to a period derivative of

[TABLE]

This is roughly $3\sigma$ lower than the value of $-12.6\pm 1.2$ ms yr*-1* found by B19. We investigated this by applying our analysis to only the transit times used by B19, which precisely reproduced their results. Therefore the difference in the period derivatives found by ourselves and by B19 is fully explained by the inclusion of additional transit times in the current study.

The quadratic ephemeris has a much lower $\chi_{\nu}^{\,2}$ than the linear ephemeris. To explore this further we calculated the Bayesian Information Criterion (Schwarz, 1978) to be 172.2 for the quadratic and 269.6 for the linear ephemeris, and the Akaike Information Criterion (Akaike, 1974) to be 165.1 for the quadratic and 264.8 for the linear ephemeris. All three statistical quantities give a significant improvement for the quadratic versus the linear ephemeris, meaning that we have detected TTVs in the WASP-4 system. The quadratic term is negative and measured to 8.4$\sigma$ significance.

The $\chi_{\nu}^{\,2}$ of our preferred quadratic ephemeris is significantly above unity. We have found this situation to occur for most transit timing studies both in our own experience (e.g. Southworth et al., 2016) and in published studies (e.g. Baluev et al., 2015; Maciejewski et al., 2018a). The 8.4$\sigma$ significance level of the quadratic term already includes an allowance for this large $\chi_{\nu}^{\,2}$. Fig. 2 shows that our new timings are consistent with published values in all cases where such a comparison can be drawn.

WASP-4 A shows starspots and these are expected to add scatter to timing measurements both on theoretical (Ballerini et al., 2012; Oshagh et al., 2013) and empirical (Maciejewski et al., 2018b) grounds. The obvious outlier at cycle $-544$ and residual $-64$ s is affected by starspots. To investigate this we rejected the six timings clearly affected by starspot crossing events (Table 1) and recalculated the ephemerides. We find that the $\chi_{\nu}^{\,2}$ drops to 1.30 for the quadratic and 1.70 for the linear ephemeris, the quadratic term is not significantly changed, and its detection significance becomes 8.2$\sigma$. We conclude that spot activity increases the scatter of the transit timings but that TTVs are still detected to high confidence.

Visual inspection of Fig. 2 suggested the possibility of a sinusoidal variation superimposed on the quadratic term. We therefore calculated a periodogram of the residuals of the best fit to the quadratic ephemeris, using the period04 package (Lenz & Breger, 2004). The periodogram contained no significant peaks, with the highest having a signal-to-noise ratio of 2.7. We find no significant peak near the period of 5.1 d tentatively identified by Baluev et al. (2015), so do not confirm this detection of sinusoidal variation. Further observations are needed to refine these conclusions.

5 Possible causes of the transit timing variation

In the previous section we have demonstrated that WASP-4 exhibits a quadratic variation in its transit times with a significance of 8.4$\sigma$. We now discuss possible causes of this effect.

5.1 Starspots

Transit times measured from light curves exhibiting spot crossing events can be biased away from the true time of midpoint because the fitted model (transit of a planet across a limb-darkened but otherwise uniform-brightness star) is no longer adequate (Ballerini et al., 2012; Oshagh et al., 2013). The bias depends on where the spot crossing event is, with those near the start or end of a transit having the strongest effect (Oshagh et al., 2013). As activity levels on stars undergo cycles lasting years to decades (e.g. Baliunas et al., 1995), these effects could conceivably mimic a slow TTV. We have visually inspected all transits presented in the current work in order to find those with clear spot crossing events. We have previously found visual methods to perform well in the identification of spot anomalies in transit light curves (Močnik et al., 2017). We have indicated the transits with spot crossing events on Fig. 2, and find no clear correlation with the residuals versus the quadratic ephemeris fit. The quadratic term in the ephemeris is also detected to almost the same significance if these timings are rejected from the ephemeris fit.

From Fig. 1 we see that starspots are preferentially detected in WASP-4 towards the start or end of totality in the transit light curve. This is because the starspots move with the stellar surface at an approximately constant angular speed, and thus at a variable linear speed when projected into the plane of the sky. The starspots move faster when they are near the projected centre of the star, so are less likely to be found there. Conversely, smaller starspots may be preferentially found away from the limb of the star because limb darkening will attenuate their photometric signal and make their detection less likely. Both effects should be accounted for in detailed statistical studies of these phenomena.

Watson & Dhillon (2004) investigated the effect of starspots on timing measurements due to the Wilson depression: the phenomenon that starspots are slightly depressed into the surface of the star. They found that this effect was sufficient to cause a change in the measured time of an eclipse at the level of a few seconds for a white dwarf plus M-dwarf binary. The effect will be less than a second for WASP-4 due to the larger radius ratio, so is not capable of causing the TTV we have found.

5.2 Applegate mechanism

Changes in the quadrupole moment of the star (Applegate, 1992) have been proposed as a possible cause of TTVs in planetary systems. Watson & Marsh (2010) calculated the amplitude of the TTVs for WASP-4 for three scenarios concerning the length of the activity cycle in the host star. The largest TTV amplitude was found to be 15.3 s for a 50 yr activity cycle, so is a factor of six too small to explain the observed TTVs.

5.3 Orbital decay

Tidal effects will cause the majority of known hot Jupiters to merge with their host stars on a timescale of several Gyr (Levrard et al., 2009; Jackson et al., 2009). This orbital decay will initially appear as a continuous decrease in their orbital period, causing a quadratic variation in the transit times. We followed the approach of Patra et al. (2017) to calculate the modified stellar tidal quality factor for the WASP-4 system, finding $Q_{\star}^{\prime}=10^{4.58\pm 0.03}$.

A wide range of possible $Q_{\star}^{\prime}$ values have been found in previous studies, both theoretical and observational. The canonical value of $Q_{\star}^{\prime}$ is $10^{6}$ (e.g. Ogilvie & Lin, 2007), a value of $10^{5.5}$ was found to be a good match to a sample of extrasolar planets (Jackson et al., 2008), and a theoretical study by Essick & Weinberg (2016) yielded $Q_{\star}^{\prime}\approx 10^{5}$ to $10^{6}$ for short-period hot Jupiters. But other studies disagree: theoretical work by Penev & Sasselov (2011) constrained $Q_{\star}^{\prime}$ to be $10^{8}$–$10^{9.5}$, an observational study by Penev et al. (2012) found $Q_{\star}^{\prime}>10^{7}$, and Collier Cameron & Jardine (2018) used the distribution of orbital separations of the known population of transiting hot Jupiters to deduce $Q_{\star}^{\prime}=10^{8.26\pm 0.14}$.

The value of $Q_{\star}^{\prime}$ is also expected to depend on the nature of the tidal perturbation and the internal structure of the star (Ogilvie, 2014; Penev et al., 2018), on orbital period (Barker, 2011), and on planet mass (Barker & Ogilvie, 2010). Barker (2011) found a dependence on orbital period of $Q_{\star}^{\prime}\propto P^{14/5}$, giving a value of approximately $10^{5.3}$ for the WASP-4 system.

A tidal quality factor small enough to cause the TTV observed in WASP-4 is smaller than any determined in the observational and theoretical studies listed above. However, it is sufficiently close to the lower envelope of $Q_{\star}^{\prime}$ values, and this envelope is sufficiently large, that it is not currently possible to rule out orbital decay as a cause of the TTV. Further work is needed to better understand the possible values of $Q_{\star}^{\prime}$ and its dependence on system parameters such as planet mass and orbital period.

5.4 Light-time effect

The presence of a third component in the system on a wider orbit than WASP-4 b will cause a change in the measured transit timings due to the varying distance of the transiting system from the observer coupled with the finite speed of light.

As our data are adequately explained by a quadratic function, we can provide only speculative limits on the properties of a sinusoid that might also provide a good fit to the timings. To do so, we fit a linear ephemeris to all timings obtained before JD 2456600 and found that the most recent data arrived $\sim$100 s earlier than predicted by this ephemeris. As lower limits we adopted a period of triple the time interval covered by all of the timings, and an amplitude of 300 s. This yields a mass function of $2\times 10^{16}$ kg, giving a minimum mass of the putative third component of \sim$$2\times 10^{-5}$$\,{\rm M}_{\odot} (0.02$\,{\rm M}_{\rm Jup}$). An object of this mass would have easily evaded spectroscopic detection, and at a minimum orbital separation of 0.6 au it would be much too close for detection via high-resolution imaging.

If present, this third object would imprint an orbital motion of amplitude $\sim$500 m s*-1* on the radial velocity (RV) of the star. Knutson et al. (2014) searched for changes in the systemic velocity of WASP-4 A, finding an insignificant value of $\dot{\gamma}=-0.0099^{+0.0052}_{-0.0054}$ m s*-1* d*-1*. The 3$\sigma$ upper limit on this is 5.9 m s*-1* yr*-1*, so would be 36 m s*-1* over the 6-year time interval covered by the RVs and 69 m s*-1* over the 11.7-year time interval covered by the transit timings. The lowest possible RV change over a 6-year time interval is 67 m s*-1*, two times larger than allowed by the 3$\sigma$ upper limit on $\dot{\gamma}$, so the third-body interpretation is moderately inconsistent with the data.

Additional possibilities that could be invoked to bring the third-body hypothesis into better agreement with the data are orbital eccentricity (which must be sufficiently low to yield a roughly quadratic variation over the 11.7-year interval covered by the available timings), an inclined orbit (which will not help much because it would require a larger third body mass to obtain the same amplitude of the light-time effect), and a longer orbital period (which would also require a more massive third body). We conclude that the light-time effect is not a good explanation for the observed TTVs. Further analysis and RV measurements of the host star would help to strengthen this conclusion.

5.5 Apsidal motion

Patra et al. (2017) attempted to explain the quadratic orbital ephemeris of WASP-12 through apsidal precession, finding that it was consistent with the data but disfavoured relative to an orbital decay model. It is difficult to apply this model to WASP-4 because it requires measurements of the midpoint of occultations as well as transits over a significant time interval. Whilst some occultation timings are available (Beerer et al., 2011; Cáceres et al., 2011; Zhou et al., 2015) they lack the precision and time sampling necessary to provide useful constraints on apsidal motion in the system (see also B19). We recommend that new high-precision observations of occultations in the WASP-4 system are obtained.

Apsidal motion can only occur in eccentric systems. Triaud et al. (2010) found $e<0.018$ for WASP-4 and Beerer et al. (2011) found $e\cos\omega=0.00030\pm 0.00086$: these do not constrain eccentricity sufficiently to rule out apsidal precession. One problem for this hypothesis, however, is the need for the system to have a non-zero eccentricity despite the presence of strong tidal effects in this short-period system containing a cool star with an extensive convective envelope.

6 Summary

We have presented 22 times of mid-transit for the WASP-4 transiting planetary system that confirm the early arrival of transits noticed in TESS data (B19) and are consistent with a quadratic ephemeris. A constant-orbital-period model is rejected at the 8.4$\sigma$ significance level, making WASP-4 b only the third hot Jupiter with TTVs detected to high significance.

We have explored possible explanations for the observed TTVs. Instrumental timing effects cannot be blamed because there is good agreement between different telescopes in the same observing seasons, and they are also not apparent in any of the tests we have performed. Biases due to spot crossing events in the light curves have a negligible effect on our results so can be ruled out as the source of the TTVs. The Applegate mechanism gives TTVs a factor of six smaller than observed. The light-time effect due to a third body struggles to match our results. Apsidal precession is plausible and deserves further investigation, in particular the observation of secondary eclipses so the orbital period found from transits can be compared to that found from occultations.

Orbital decay is a plausible but not favourable option for causing the observed TTV. The tidal quality factor needed, $Q_{\star}^{\prime}=10^{4.58\pm 0.03}$, is smaller than those predicted from theoretical arguments or measured from the population of known transiting planets. However, the envelope of possible values of $Q_{\star}^{\prime}$ is broad and close to the value found for WASP-4, so it is not currently possible to discount this hypothesis.

A program of transiting timing measurements of a wide variety of planetary systems would help by providing observational constraints on orbital decay as a function of stellar type, planet mass, and orbital period. As $Q_{\star}^{\prime}$ is expected to depend on the system parameters, it also would be interesting to perform TTV analyses on other systems similar to WASP-4. We used the TEPCat333TEPCat (the Transiting Extrasolar Planet Catalogue) is available at http://www.astro.keele.ac.uk/jkt/tepcat/ catalogue (Southworth, 2011) to identify HATS-2 and WASP-64 as the best candidates. These were found 5–6 years later than WASP-4 so their transit timings cover a smaller time interval, but this could be compensated for by obtaining a larger number of high-quality transit observations in the near future.

We also advocate obtaining new observations of the WASP-4 system, and will be keeping it on the target list of the MiNDSTEp transit monitoring program. New transit times would refine and extend the baseline of the quadratic ephemeris. Further timings of secondary eclipses would allow the apsidal motion hypothesis to be ruled in or out. RV measurements of the planet host star would give stricter limits on the change in systemic velocity of the system. WASP-4 and WASP-12 stand out as the only hot Jupiters for which unexplained TTVs have been detected: as this phenomenon is uncommon it is reasonable to propose unusual or unlikely scenarios for such systems.

Acknowledgements

This paper uses data from the South African Astronomical Observatory (SAAO) under programmes Southworth-2016-05-1.0-m and Southworth-2019-01-40-inch-228, and from ESO under programme ID 088.C-0204 (P.I. Tregloan-Reed). The reduced light curves presented in this work will be made available at the CDS (http://cdsweb.u-strasbg.fr/) and at http://www.astro.keele.ac.uk/jkt/. We are grateful to Luke Bouma, Adrian Barker and an anonymous referee for very useful comments that helped to significantly improve our paper. This project was partially funded by the MINEDUC-UA ESR project ANT 1795. The following internet-based resources were used in research for this paper: the ESO Digitized Sky Survey; the NASA Astrophysics Data System; the SIMBAD database operated at CDS, Strasbourg, France; and the ar$\chi$iv scientific paper preprint service operated by Cornell University.