© ESO 2021

1. Introduction

Planetary nebulae (PNe) are the product of low- and intermediate-mass stars that underwent extensive mass loss during the asymptotic giant branch (AGB) phase. The simple explanation of PNe is that they are the product of stars between 1 and 8 solar masses that are nearing the ends of their lives. These stars shed their outer layers, growing brighter and hotter, and ionising the newly formed nebula before ultimately cooling and fading into white dwarfs. While this explanation may indeed hold for some PNe, it has become increasingly clear in recent decades that many of these stars do not die alone: A significant fraction, if not the majority, of PNe are produced by systems that have undergone binary interactions (Miszalski et al. 2009a; Boffin & Jones 2019). Particularly compelling evidence comes from the diverse range of morphologies observed in PNe (Balick & Frank 2002), the most dramatic of which are thought to come from the presence of a close binary companion (Miszalski et al. 2009b; Jones & Boffin 2017). Some authors have even suggested that binarity might be a requirement for the formation of a PN, which would resolve the tension between population synthesis models and the unexpectedly low number of PNe observed in the Milky Way (Moe & De Marco 2006).

Close binary central stars of PNe are not only relevant for the study of PNe: They also place important constraints on common-envelope (CE) theory (Iben & Livio 1993; Yungelson et al. 1993), a poorly understood evolutionary phase that is critical for the formation of a wide variety of astrophysical phenomena, including cataclysmic variables, low-mass x-ray binaries, type Ia supernovae, mergers, and gravitational wave sources (Jones 2020). A statistically significant sample of close binary systems is crucial for observationally constraining the orbital period distribution, which, in turn, is an important probe of the CE efficiency (Toonen & Nelemans 2013).

Such a sample is lacking: The identification of close binary central stars has been a painstaking process, with only 15 such systems revealed in the 30 or so years following the discovery of the first example (Boffin & Jones 2019, for a historical overview) and the current total standing around 70. Binarity manifests itself in various ways. Widely separated pairs can be detected visually, particularly from space (Ciardullo et al. 1999; Benetti et al. 2003; Liebert et al. 2013), and bright but close companions produce blended spectra (which can sometimes be confirmed by long-term, high-resolution radial velocity monitoring; Jones et al. 2017). Radial velocity variations are the most conclusive evidence, but are often difficult to measure due to contamination from bright nebula lines and spurious variability stemming from post-AGB winds (De Marco et al. 2004). The most dramatic photometric variability comes from eclipses (requiring a lucky alignment of the orbital plane, e.g., Jones et al. 2019) and/or the irradiation of a cool companion (e.g., Munday et al. 2020). In very close systems, tidal distortion can result in the presence of ellipsoidal modulation in the observed light curve (e.g., Santander-García et al. 2015). Detection of photometric variations is also subverted by the presence of a bright nebula and aggravated by variable seeing conditions for ground-based monitoring (Jones et al. 2015). Space-based photometric monitoring has shown promising results, with variability in PN central stars having recently been identified in Kepler/K2 (De Marco et al. 2015; Jacoby et al. 2020) and TESS (Aller et al. 2020), though the high sensitivity of these missions is somewhat offset by the low spatial resolution, which dilutes the signal, particularly in the presence of nearby stars or bright nebulosity.

The European Space Agency’s Gaia mission (Gaia Collaboration 2016) offers great promise as it tracks not only positional but also photometric variations, and does so for the whole sky, including most known PN central stars (Chornay & Walton 2020). Gaia offers high spatial resolution and long time baselines at the cost of irregular sampling and a long wait for epoch photometry to be available for the complete set of sources. This paper illustrates the potential of the Gaia mission for the discovery and characterisation of new binary PN central stars.

2. Variability from mean photometry

Over the course of Gaia’s mission, any given source is observed many times. The pattern of the observations is dictated by Gaia’s scanning law (Gaia Collaboration 2016), and the frequency at which a given source is observed depends on its position in the sky. The scanning law is such that a source should be observed at least six separate times over the course of a year, though due to Gaia’s two telescopes and its rate of rotation, each of those so-called visibility periods generally has at least two observations of the source, separated by a few hours (Eyer et al. 2017).

A crossing of a source across Gaia’s focal plane is called a field-of-view (FoV) transit. During an FoV transit a source’s flux is measured multiple times across different charge-coupled devices (CCDs) and through different passbands. Gaia has three passbands: the main, wide G band, and the blue G_BP and red G_RP bands. In the G band there are typically eight or nine CCD transits per FoV transit; in G_BP and G_RP there is one at most.

For most sources, Gaia’s intermediate data releases include only the mean photometry in each band, that is, the mean flux $\overline{I_{\mathrm{s}}}$ and its uncertainty ${\sigma }_{\overline{I_{\mathrm{s}}}}$. In the calculation of these quantities individual CCD flux measurements are weighted by their nominal precisions (Eqs. (6) and (7) in Carrasco et al. 2016). However, the uncertainty in the mean is based on the observed measurement scatter, so it may not match that expected from the uncertainties of the individual fluxes. In particular, if a source is variable, this will manifest as increased scatter and a higher uncertainty. If variability is the dominant source of the scatter, the amplitude of that variability can be derived from the published uncertainty.

2.1. Variability amplitude

The uncertainty in the mean flux of a source published in the Gaia catalogue is the weighted standard deviation of the individual flux measurements divided by the square root of the number of CCD transits N_CCD (corresponding to the phot_g_n_obs column). The weights are not published, but if we assume that the weights are equal, we can recover the standard deviation of the fluxes in the G band as $\sqrt{N_{\mathrm{CCD}}}{\sigma }_{\overline{I_G}}$, where now $\overline{I_G}$ and ${\sigma }_{\overline{I_G}}$ are the mean and uncertainty for the G -band flux.

Following Mowlavi et al. (2021), Belokurov et al. (2017), and others, we adopt the variability statistic

$$\begin{array}{c}\hfill A_{\mathrm{var},G}\equiv \sqrt{N_{\mathrm{CCD}}}\frac{{\sigma }_{\overline{I_G}}}{\overline{I_G}},\end{array}$$(1)

which is simply the standard deviation of the individual CCD fluxes as a proportion of the mean flux. This is related to the scatter in the magnitude as follows. We recall that the Gaia G -band magnitude $G=-2.5{log}_{10}\overline{I_G}+G_0$, where G₀ is the zeropoint. For a scatter in the individual flux measurements that is small relative to the mean flux (A_var, G ≪ 1), the standard deviation σ_m of the corresponding magnitudes is

$$\begin{array}{c}\hfill {\sigma }_m=\frac{2.5}{ln10}{A}_{\mathrm{var},G},\end{array}$$(2)

using the properties of the derivative (i.e. the first-order Taylor polynomial).

An even more convenient relation comes from assuming that the variation is sinusoidal, as it often approximately is for close binary systems. In this case, ${\sigma }_m=A/\sqrt{2}$, where A is the sinusoid semi-amplitude. Equivalently, the range 2A is then

$$\begin{array}{c}\hfill \mathrm{range}G=2\sqrt{2}{\sigma }_m\approx 3.07A_{\mathrm{var},G},\end{array}$$(3)

where σ_m is as in Eq. (2). For the full derivation and a deeper discussion, see Mowlavi et al. (2021).

2.2. Variability significance

The intrinsic measurement uncertainty varies significantly with magnitude, increasing for fainter sources due to photon statistics as well as including noticeable features at other magnitudes related to differences in data acquisition and calibration. It is important to ascertain whether the apparent variability that we obtain is actually significant. If we had the per-FoV transit uncertainties, we could determine an excess noise statistic analogous to the unit weight error for astrometric measurements. However, we do not have those, and, moreover, the errors are heteroscedastic, depending not only on source colour and magnitude, but varying temporally and spatially due to factors such as straylight conditions and sky background. Nonetheless, the considerations are broadly similar to those discussed in Lindegren (2018) in the context of Gaia’s re-normalised unit-weight error (RUWE) statistic for astrometry.

With the aim of obtaining a crude significance statistic, we assumed that there is an intrinsic uncertainty in the flux measurement for a single FoV transit that depends on a source’s magnitude and colour, σ_G = σ_G(G, C) (where colour is the magnitude difference between the G_BP and G_RP bands, G_BP − G_RP). We refer to the standard deviation of the observed fluxes across different FoV transits for a given source using σ_G, obs; in general, this will only match σ_G in the limit of a large number of samples.

The arrangement of CCDs on Gaia’s focal plane results in sources having an average of 8.86 G-band flux measurements per FoV transit (Jordi et al. 2010). As CCD transits occur sequentially and the FoV transit time is tens of seconds (Riello et al. 2018), intrinsic source variability on any much larger timescale contributes to flux scatter only across different FoV transits. Thus we focus on FoV rather than CCD transits, and define the effective number of samples k ≡ N_FOV ≡ N_CCD/8.86.

For constant sources, values of ${\sigma }_{G,\text{obs}}^2$ will follow a χ² distribution with k − 1 degrees of freedom,

$$\begin{array}{c}\hfill k{\sigma }_{G,\mathrm{obs}}^2\sim {\sigma }_G^2{\chi }_{k-1}^2.\end{array}$$(4)

This will not match the distribution that we observe, not only because of the simplifying assumptions, but because all stars are variable at some level. However, the mode of the distribution should not be too affected by the addition of these intrinsically variable sources, so, as the mode of ${\chi }_{k-1}^2$ is (k − 1)−2, we can find σ_G as

$$\begin{array}{c}\hfill {\sigma }_G^2=\mathrm{mode}\left(\frac{k{\sigma }_{G,\mathrm{obs}}^2}{(k-1)-2}\right),\end{array}$$(5)

assuming that the empirical distribution of ${\sigma }_{G,\text{obs}}^2$ over which the mode is taken is appropriately smoothed.

For sufficiently large k, the ${\chi }_{k-1}^2$ distribution can be approximated as a normal distribution with matching mean k − 1 and variance 2(k − 1),

$$\begin{array}{c}\hfill \frac{k{\sigma }_{G,\mathrm{obs}}^2}{{\sigma }_G^2}\sim {\chi }_{k-1}^2\approx \mathcal{N}(k-1,2(k-1)).\end{array}$$(6)

This can be related to a standard normal by subtracting the mean and dividing by the square root of the variance,

$$\begin{array}{c}\hfill \sqrt{\frac{k-1}{2}}(\frac{k{\sigma }_{G,\mathrm{obs}}^2}{(k-1){\sigma }_G^2}-1)\sim \mathcal{N}(0,1).\end{array}$$(7)

The term on the left is a dimensionless measure of significance, which can be further approximated for convenience as a natural logarithm using the relation x ≈ ln(1 + x) for small x (as the part inside the parenthesis is close to zero). Furthermore, an additional magnitude- and colour-dependent scale factor can be added to ensure that the variance is indeed 1 (in practice, to ensure that the left side of the distribution has variance 1, as the right side includes most variable sources). This gives

$$\begin{array}{c}\hfill D\equiv \frac{1}{D_0}\sqrt{\frac{k-1}{2}}ln\left(\frac{k{\sigma }_{G,\mathrm{obs}}^2}{(k-1){\sigma }_G^2}\right),\end{array}$$(8)

where D is the variability significance and D₀ = D₀(G, C) is the scale factor. It is straightforward to compare this form of the expression to the logarithm of our variability amplitude A_var, G.

In Eq. (8) the significance of any observed variability increases with the number of FoV transits k. However, as k increases, the simplifying assumptions (rather than statistical properties) will increasingly contribute to the observed scatter.

We determined σ_G and D₀ empirically using a representative sample of Gaia sources. We binned sources by magnitude and colour and calculated σ_G and D₀ using the distribution of flux uncertainties in each bin. Interpolating between these values yields two-dimensional functions of magnitude and colour that we used to compute the variability significance for individual sources.

3. Results

We applied our methodology to the CSPN catalogue of Chornay & Walton (2020). This catalogue was generated using an automated method to match CSPNe with sources in Gaia data release 2 (DR2), and includes a reliability parameter (between 0 and 1) that indicates the certainty of the match. We expect that Gaia’s early data release 3 (EDR3) will not significantly change the results; in particular, it does not contain any new epoch photometry or variability analysis that can be used to validate our assumptions, and any new detections will be for faint sources where only the most extreme variability would be detectable. We join our catalogue with a compilation of known and suspected binary CSPNe from literature¹, excluding one extragalactic PN and six objects identified as non-PNe in the Hong Kong/AAO/Strasbourg Hα (HASH) PN catalogue of Parker et al. (2016), while adding PHR J1040-5417 (Hillwig et al. 2017).

3.1. Known binaries

Figure 1 shows the variability amplitude plotted against Gaia G magnitude for the sources in our sample, colour-coded by binary type. The sample has been filtered by CSPN match certainty (reliability > 0.5) and PN angular size (radius ≥10″; see Sect. 3.2 for motivation), giving 683 objects in total (see Table 1).

Fig. 1. Variability amplitude vs. Gaia G magnitude for CSPNe in our sample. Filled markers are known or suspected binary systems of different observational types; see text for description. Lines indicate constant variability significance D for different numbers of FoV transits N.

Non-variable binaries (blue) include those with composite spectra and those detected (or hypothesised) based on infrared excess. Visual binaries (green) have typically been discovered through high resolution Hubble Space Telescope imagery (and appear blended in ground-based imagery), but do not have measured radial velocity or photometric variation. Long-period binaries (orange) have composite spectra and measured periods (using radial velocity variations) on the order of years. Some exhibit shorter timescale variability due to rotation. Short-period or close binaries (red) generally exhibit the greatest photometric variability. These have typical periods shorter than one day, with only a few between 2 and 20 days. Half-filled red circles are possible short-period binaries whose nature is still uncertain: the recent discoveries in TESS from Aller et al. (2020), objects designated as in need of further study by Miszalski et al. (2011), and Pa 5 from Kepler (De Marco et al. 2015). A few short-period binaries had also been noted previously as visual binaries or composite spectra; in these cases, the short-period binarity takes precedence. The remainder of the sample (objects of unknown binarity) are shown as empty circles. Markers outlined in black are objects that have low astrometric error (RUWE < 3; see Sect. 3.2) and significant variability (D > 8; see Sect. 3.3).

The lines in the plot are lines of constant variability significance D = 0 and D = 8 (see Eq. (8)) for a source with a representative CSPN colour, that is, G_BP − G_RP = 0 (this value, which is the 44th percentile of the colours of the CSPNe in the plot, is also used for sources with no colour information). The D = 8 line is split into two different lines for 15 and 30 FoV transits, which are typical counts representing the 21st and 63rd percentile of N_FOV for CSPNe in the sample. D scales with the square root of N_FOV (k in Eq. (8) and preceding equations). The non-monotonic pattern in typical photometric uncertainty is evident, and is similar in shape to the magnitude error distribution in Fig. 9 of Evans et al. (2018), with peaks around G = 11 and 13 due to different instrumental configurations (gates and window classes; Carrasco et al. 2016) affecting calibration. Towards fainter magnitudes, where photon noise dominates, the errors increase more predictably. Though not shown in the plot, the intrinsic scatter for sources with extreme colours (e.g., very blue CSPNe) also tends to be higher, particularly at brighter magnitudes, though the overall shape of the magnitude dependence is the same.

In Fig. 2 we further break down known close binary systems by the type of variability. The criteria for inclusion are loosened relative to Fig. 1 in order to increase the sample size (see counts in Table 1). The greatest variability is evident in irradiated binaries. Ellipsoidal modulation generally produces a smaller effect, while most binaries that have been detected solely through radial velocity variation (spectroscopic binaries, shown as upward pointing triangles) have no significant photometric variability in Gaia, being for instance double-degenerate systems such as Fg 1 (Boffin et al. 2012) and NGC 2392 (Miszalski et al. 2019). Of the systems with low Gaia variability, only AMU 1 has been previously identified as a photometric variable, with Kepler photometry revealing low-amplitude modulation due to relativistic Doppler boosting (De Marco et al. 2015).

Fig. 2. Same as Fig. 1, but for short-period binaries only and with no minimum radius cut. The relaxation of the radius requirement means that not every point on this plot has a corresponding red circle in Fig. 1. Different marker shapes indicate different means of binarity detection, with filled markers being eclipsing (possibly in addition to another type of photometric variation).

These figures demonstrate the power of the scatter in the Gaia mean photometry. Comparing Figs. 1 and 2, we can see that all short-period binaries in the sample of extended PN with CSPNe brighter than G = 19 that have ground-based detections of photometric variability (i.e. that are not solely spectroscopic binaries) also have significant variability in Gaia, with variability amplitudes corresponding to peak-to-peak amplitudes of at least 0.1 mag in the G band (around the limit of what is detectable from the ground). The only confidently known extended binary fainter than G = 19, MPA J1759-3007, is visibly separated from the trend, but has a low number of FoV transits and narrowly fails to meet the D > 8 significance cut.

3.2. Other sources of photometric scatter

Without time-series photometry it is difficult to distinguish between genuine variability (possibly related to binarity) and unusually high photometric scatter. Being bright, extended objects, PNe are particularly likely to cause problems for Gaia. In this section we consider possible sources of spurious variability.

Only single stars are modelled in the Gaia DR2 astrometric solution. The quality of that fit is indicated by the RUWE statistic (Lindegren 2018). Unresolved binary systems can introduce wobble in the location of the flux maximum, which will contribute to a poor fit and a high value of RUWE (Belokurov et al. 2020), though this is most noticeable for sources near us (within a few hundred parsecs). For visible binaries separated by less than 0.4″, astrometric issues can also be caused by the two components of a binary being detected separately and later incorrectly merged into a single source (Antoja et al. 2018), in which case the duplicated_source flag will be set. This merging is also likely to manifest as spurious photometric variability if the components have different magnitudes. Partially resolved binary systems could also cause issues for Gaia’s point or line spread function fitting, introducing a scan angle dependence to the flux and position measurements.

The top panel of Fig. 3 shows RUWE versus the variability amplitude significance, with the same colouring scheme as in Fig. 1. All known extended short-period binary systems (red points) have RUWE < 3, with most having RUWE < 1.4, 1.4 being the rather conservative cutoff for well-behaved sources suggested by Lindegren (2018). These values suggest that the photometric scatter is real rather than due to cross-matching or fit issues. Several visual binaries (green) have significant variability and high RUWE values. These pairs are typically separated by around 0.3″, a separation at which they can be resolved by Gaia, but then are merged into a single source. This merging is likely responsible for both the high photometric scatter and poor astrometric fit.

Fig. 3. Quantities potentially associated with spurious variability in Gaia photometry. Markers are different binary types as in Fig. 1. Vertical lines correspond to the lines of constant variability significance in Fig. 1. The significance of the horizontal lines is described in the text. Duplicated sources in Gaia are indicated in the upper plot with pluses. PNe in the upper and middle plots are limited to a minimum radius of 10″ (corresponding to the selection in Fig. 1), while the lower plot has no such restriction.

Several visual binaries also have relatively high photometric excess factors (middle panel of Fig. 3), as the window used for G_BP and G_RP photometry captures both members of the pair. However, the distinction between well-behaved and poorly behaved sources is less clear than for RUWE, as other features such as crowded fields and bright nebulae can also contaminate the G_BP and G_RP photometry without necessarily introducing spurious variability in the G band. The suggested excess factor threshold from Evans et al. (2018) for an object with G_BP − G_RP = 0 is shown for reference, though it is too conservative for our application, and we did not use any cut based on the excess factor.

The lower panel of Fig. 3 shows the relation between variability significance and nebula radius. There is a notable excess of high-variability objects at smaller radii. This is not surprising, as these nebulae tend be brighter, and thus introduce an additional source of sky background. Compact nebulae can also cause the sorts of astrometric issues described previously. Known binaries at smaller radii do exhibit significant variability, but precisely selecting them is difficult. Thus we adopted a minimum radius of 10″ in our selection.

3.3. New candidate binaries

Markers outlined in black in Fig. 1 are objects that have significant Gaia variability (D > 8) and low astrometric error (RUWE < 3), in addition to meeting reliability > 0.5 and radius ≥10″ criteria common to all objects in the extended sample. We chose the criteria, in particular the threshold D = 8, as a balance between precision and recall, with the aim of selecting known short-period binaries (red points in Fig. 3). The 92 objects meeting these criteria represent 13% of that sample overall, but include 27 short-period binaries (77% of those in the sample), as well as 4 long-period binaries exhibiting shorter timescale variability, and one visual binary with a borderline RUWE value.

This leaves 58 objects that are not previously known or suspected binaries but exhibit significant variability in Gaia. Of these, 43 have log variability amplitude >  − 1.5, suggestive of a level of variability that, if real, should be recoverable from ground-based follow-up observations. Confirming and measuring the periods of these systems, even at a moderate success rate, would dramatically increase the sample of known binary CSPNe. The list of these objects is given in Table A.1.

Such a wealth of previously undiscovered binaries is not unexpected. An unbiased ground-based survey, OGLE-III, previously found a close binary fraction of 12–21% (Miszalski et al. 2009a). Our sample of extended PNe (radii ≥10″) with well-behaved astrometry (RUWE < 3) and bright central stars (G < 18) contains 318 PNe. Of these, 9–11% are known short-period binaries. If we were to assume that all of the detected variables with log variability amplitude ≥1.5 are indeed binaries (i.e. trying to recreate the limits of the OGLE survey), we would find a binary fraction of 18%, consistent with the range from OGLE. The estimate depends on the particular thresholds chosen and misses spectroscopic binaries, but this consistency supports the hypothesis that many of the newly discovered variables are indeed binaries.

3.4. Nebular chemistry

One strong indicator of PN central star binarity is a high abundance discrepancy factor (ADF; Wesson et al. 2018). This is the ratio of abundances of the same species derived from recombination lines of a photoionised nebula to those derived from collisionally excited lines. It has a median value of 2.6 for PNe. Extreme values of the ADF (ADF > 5) in PNe are associated with close binary central stars with a period of the order of a day or less. Showing that some of the variables in our sample exhibit such extreme ADFs would strengthen the link.

We joined our catalogue with a compilation of literature ADF values² and again required reliability > 0.5, but relaxed the angular size cut to include more compact PNe (radii ≥5″). The 50 objects that passed these cuts (out of the 123 Galactic PNe with literature ADFs) are shown in Fig. 4. There are 21 PNe exhibiting extreme ADFs, of which 12 are already known or suspected to harbour close binaries.

Fig. 4. Extreme ADF objects (indicated with crosses) on the variability vs. magnitude plot. Labelled objects are PNe with extreme ADFs that are not previously known binaries. Marker colours indicate type of binarity if known, and the grey lines are the same lines of constant variability significance shown in Fig. 1. Likewise, markers outlined in black have the same meaning as in Fig. 1.

As in Fig. 1, the lines indicate constant variability significance D = 0 and D = 8, and circles with black outlines are sources with both RUWE < 3 and D > 8. The sources above the D = 8 lines without black outlines have high RUWE values; the bias of the ADF sample towards compact, bright PNe, which are more amenable to spectroscopy, likely contributes to the larger proportion of astrometric issues relative to the sample in Fig. 1. These issues could in turn be responsible for some of the observed photometric scatter (as it seems to be for some sources in the top panel of Fig. 3), though the observed scatter still serves as an upper bound.

With this in mind, it is reassuring that of the extreme ADF sample, only NGC 40 does not show at least some variability. Along with NGC 40, NGC 2022 and M 1-42 have so far eluded ground-based detections of periodic variability. From Fig. 4 we can see that variation in the CSPN of NGC 2022 seems likely to be genuine (D = 5.7, RUWE = 1.4) but on the threshold of ground-based detectability, while for M 1-42 with high variability significance D = 13.3 but RUWE = 3.0, the situation is less clear. Of the remaining extreme ADF objects, four have significant variability and low RUWE (though most of their amplitudes are again likely to be difficult to detect from the ground), while NGC 2440 has an extreme RUWE value that makes its apparent variability inconclusive, and H 1-41 has a moderate variability significance of D = 6.8, which may be improved in the future with more observations.

Thus, from the Gaia data alone, we can conclude that NGC 40 is unlikely to have strong photometric variability (though that of course does not preclude it from being a binary), but only state for other extreme ADF objects that their photometric scatters are consistent with variability caused by binarity. Light curves are needed to ascertain whether there is a periodic component to the scatter of these objects, and whether the period matches the expectation from Wesson et al. (2018), who found that extreme-ADF PNe have central stars with periods shorter than 1.15 days.

4. Light curves from epoch photometry

Characterisation of variable sources in Gaia DR2 is limited to about 550 000 sources that were classified as VARIABLE by Gaia’s own analysis (Holl et al. 2018). Epoch photometry (flux measurements for a source at each FoV transit) is published for all of these sources, as well as a table (vari_time_series_statistics) of summary statistics derived from that photometry. About 3000 sources exhibiting short timescale variability additionally have frequency information published in the vari_short_timescale table.

However, completeness is low: Only four of the known short-period binary CSPNe are classified as VARIABLE, even though many more are known to exhibit strong photometric variability, and only two of them have frequencies recovered by Gaia. Additionally, four other sources are classified as VARIABLE that are not previously known binaries. All eight of these VARIABLE sources were also identified in our selection as significantly variable based on their flux uncertainties. Though they represent only a fraction of such sources, their epoch photometry offers a window through which to probe the nature of their variability, as well as the assumptions that went into our analysis.

4.1. Frequency extraction

We performed a Lomb-Scargle analysis (using the implementation of VanderPlas & Ivezić 2015) on the eight G-band light curves in order to extract periodic variability, with a grid of 10⁶ frequencies logarithmically spaced between 10^−1.5 and 10^1.5 days⁻¹. We selected the frequency at the peak with the highest power, and fitted the light curves to a sinusoid

$$\begin{array}{c}\hfill G(t)=Asin(2\pi ft+\varphi )+c,\end{array}$$(9)

where A and c are the amplitude and offset in mags, f is the frequency in days⁻¹, t is the barycentric Julian date (JD) in barycentric coordinate time (TCB) relative to 2010-01-01T00:00:00, and ϕ is the phase. We determined these parameters through least-squares optimisation, allowing the frequency to vary within a small range close to the highest power frequency extracted from the periodogram.

4.2. Results

Table 2 contains the parameters of the best fits. Derived phase-folded light curves and residuals are shown in Fig. 5 for known binaries and Fig. 6 for new candidates. All eight VARIABLE sources have good fits, and three out of the four derived periods for known binary systems also match literature values.

Fig. 5. Epoch photometry (left), derived periodograms (left centre), phase-folded light curves (right centre), and residuals (right) for previously known binaries classified as VARIABLE by Gaia. Red ticks on the periodograms indicate the highest peak, and horizontal dashed lines are false-alarm levels at (from low to high) 10%, 1%, and 0.1%. Powers are normalised relative to the residuals of a constant model and lie between 0 (no improvement over constant model) and 1 (periodic model has zero residuals). Error bars are 2σ errors for visibility. An additional half-phase is shown for the light curves and residuals. The light curve and residuals for Abell 41 are folded at double the derived period (see text).

Fig. 6. Same as Fig. 5, but for new candidate binaries.

One exception is Abell 41, where our derived period is half that in literature (Jones et al. 2010). The variability of Abell 41 is believed to be due to ellipsoidal modulation and partial eclipses, so it has two minima per orbit, and the shape of its folded light curve is different from the others in the sample and less well described by a sinusoid. We have plotted the folded light curve for Abell 41 at double our derived period to reflect this, though the difference in the depths of the minima for example is difficult to see with Gaia’s sparse sampling (cf. Fig. 2 of Bruch et al. 2001).

The known irradiated binary PHR J1040-5417 has deep primary eclipses (Hillwig et al. 2017) that are not visible in the epoch photometry; instead there is a gap at the light-curve minimum. The eclipse is deep enough that it would be close to the detection limit of Gaia DR2; however, comparison with the Gaia scanning law shows that this gap also exists in the sampling.

The four newly discovered variable CSPNe are those of Cr 1, KTC 1, Kn 124, and Pa 48. All of these PNe are relatively recent discoveries; only Cr 1 and KTC 1 (Jacoby et al. 2010) have spectroscopic confirmation of their PN status. Kn 124 and Pa 48 (Kronberger et al. 2016) do not and are only listed as possible and likely PNe in the HASH PN catalogue. None have published narrow-band imagery or morphological information.

The G -band light curves of these newly discovered variables are well described by sinusoidal fits. The shape, amplitude, and periods of the curves are consistent with irradiation effects, in which one side of a cool companion is heated by the CSPN and contributes to the flux of the system most when it is directly behind the CSPN; that is, when the hot side is facing the observer. With irradiation effects there is one flux maximum per orbit, so the period of the variation is the period of the binary system. None of the systems appear to be eclipsing, and the uncertainty in G_BP and G_RP is too large to measure any colour effects. We conclude that they are likely binaries, though constraining the parameters of the systems is not possible without additional data and is beyond the scope of this work.

4.3. Relation to mean photometry

The additional data in the epoch photometry tables help to validate some of our assumptions in Sect. 2, but with a very small sample. For the eight VARIABLE sources in the vari_time_series_statistics table, the mean of phot_g_n_obs / num_selected_g_fov (i.e. N_CCD/N_FOV) is 8.89 ± 0.29, consistent with our choice of 8.86. The mean ratio of G range to A_var, G is 3.15 ± 0.10, slightly higher than but consistent with the value in Eq. (3) (though Eq. (3) applies only in the limit of uniform dense sampling).

Subtracting the best-fit sinusoids removes most of the periodic variability, though HFG 1 in particular shows some in its residuals (Fig. 5, right panels) that would be easily accounted for by adding a second higher-frequency term to its fit. The amplitude of the remaining variability (derived from the weighted variance of the residuals of the sinusoidal fits) is shown in Fig. 7. The significance of most of the residuals is low. The low residual variability of Cr 1 and KTC 1 compared to typical CSPNe is due to their relatively red central stars, which have a lower intrinsic measurement uncertainty.

Fig. 7. Initial variability (points) and residual variability (arrow tips) of Gaia VARIABLE sources after subtracting sinusoidal fits from the epoch photometry. The black lines are the same lines of constant variability significance shown in Fig. 1.

5. Conclusions

Gaia shows tremendous promise for the study of CSPNe. Though Gaia’s own variability processing is quite limited, combining the scatter in the photometry with the notion of variability significance introduced in this work allowed us to retrieve a large fraction of the binary systems that have been discovered from painstaking ground-based observations over the past decades. The newly discovered likely variables offer a promising set of targets for ground-based photometric follow-ups, building towards a statistically significant sample of post-AGB, post-CE stars. Many of their nebulae have chemistry or morphology typically associated with binarity. Selecting targets on variability rather than expectations of binarity based on other features will help us study how these features relate to binarity in a less biased manner, in addition to increasing the expected success rate.

With available Gaia epoch photometry we have shown that we can recover periods for four known binary CSPNe, and confirmed that four of the previously unknown variable sources have light curves consistent with binarity. This lends weight to our hypothesis that many of the variables that we have identified from their mean photometry are indeed binaries.

Completeness and accuracy will improve with future Gaia data releases, with more transits increasing the statistical significance of any observed excess variability. Improved photometric processing and handling of non-single systems will reduce spurious sources of photometric scatter, and ultimately, epoch photometry will enable us to confirm new binary candidates from the Gaia data alone.

Acknowledgments

We thank the anonymous referee for their comments, which have helped improve this paper. We also thank Todd Hillwig for providing a reference light curve for PHR J1040-5417. This research has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia) processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has also made use of the HASH PN database (http://hashpn.space) and of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). N. C. is supported through the Cancer Research UK grant A24042. D. J. acknowledges support from the State Research Agency (AEI) of the Spanish Ministry of Science, Innovation and Universities (MCIU) and the European Regional Development Fund (FEDER) under grant AYA2017-83383-P. D. J. also acknowledges support under grant P/308614 financed by funds transferred from the Spanish Ministry of Science, Innovation and Universities, charged to the General State Budgets and with funds transferred from the General Budgets of the Autonomous Community of the Canary Islands by the Ministry of Economy, Industry, Trade and Knowledge. D. J. also acknowledges support from the Erasmus+ programme of the European Union under Grant number 2020-1-CZ01-KA203-078200.

References

Appendix A: Table of new candidate binaries

All Tables

All Figures

	Fig. 1. Variability amplitude vs. Gaia G magnitude for CSPNe in our sample. Filled markers are known or suspected binary systems of different observational types; see text for description. Lines indicate constant variability significance D for different numbers of FoV transits N.
In the text

	Fig. 2. Same as Fig. 1, but for short-period binaries only and with no minimum radius cut. The relaxation of the radius requirement means that not every point on this plot has a corresponding red circle in Fig. 1. Different marker shapes indicate different means of binarity detection, with filled markers being eclipsing (possibly in addition to another type of photometric variation).
In the text

Fig. 3. Quantities potentially associated with spurious variability in Gaia photometry. Markers are different binary types as in Fig. 1. Vertical lines correspond to the lines of constant variability significance in Fig. 1. The significance of the horizontal lines is described in the text. Duplicated sources in Gaia are indicated in the upper plot with pluses. PNe in the upper and middle plots are limited to a minimum radius of 10″ (corresponding to the selection in Fig. 1), while the lower plot has no such restriction.

In the text

	Fig. 4. Extreme ADF objects (indicated with crosses) on the variability vs. magnitude plot. Labelled objects are PNe with extreme ADFs that are not previously known binaries. Marker colours indicate type of binarity if known, and the grey lines are the same lines of constant variability significance shown in Fig. 1. Likewise, markers outlined in black have the same meaning as in Fig. 1.
In the text

Fig. 5. Epoch photometry (left), derived periodograms (left centre), phase-folded light curves (right centre), and residuals (right) for previously known binaries classified as VARIABLE by Gaia. Red ticks on the periodograms indicate the highest peak, and horizontal dashed lines are false-alarm levels at (from low to high) 10%, 1%, and 0.1%. Powers are normalised relative to the residuals of a constant model and lie between 0 (no improvement over constant model) and 1 (periodic model has zero residuals). Error bars are 2σ errors for visibility. An additional half-phase is shown for the light curves and residuals. The light curve and residuals for Abell 41 are folded at double the derived period (see text).

In the text

	Fig. 6. Same as Fig. 5, but for new candidate binaries.
In the text

	Fig. 7. Initial variability (points) and residual variability (arrow tips) of Gaia VARIABLE sources after subtracting sinusoidal fits from the epoch photometry. The black lines are the same lines of constant variability significance shown in Fig. 1.
In the text