© The Authors 2024

1 Introduction

Planetary nebulae (PNe) are one of the most fascinating astronomical objects in the Universe. They are formed when low-to- intermediate mass stars (0.8–8 M_⊙) reach the final stage of their lives, ejecting their outer envelopes into the interstellar medium. As a result, a great diversity of PNe are observed in our Galaxy, and there are no two PNe that are exactly alike. But one thing seems clear: spherical PNe are almost the exception and they represent only one-fifth of all the PNe discovered in the Milky Way. On the contrary, the remaining 80% show a wide range of complex morphologies (Parker et al. 2006).

Explaining the precise processes that lead to the formation of these kinds of structures is a difficult task, although some progress has been made in the last decades with central star binarity emerging as the most likely culprit (see Boffin & Jones 2019, and references therein). Binary central stars of PNe (bCSPNe) appear to be key not only to explain the variety of shapes in PNe but also to respond to many other open questions in stellar evolution, such as the mass-loss rate and history, the little-understood common envelope phase, or the estimation of the binary population in general.

Thus, many of the efforts in the PN community today are devoted to the search for, and characterisation of, new binary central stars in the heart of PNe (with both short and long orbital periods), with the aim of shedding light on the formation of the striking morphologies they present. Dedicated ground-based observations of CSPNe are very costly in terms of time, so very good candidates among the large sample of PNe (more than 3000 in our Galaxy) are required to guide targeted observational campaigns. Within this context, large photometric surveys such as OGLE (Optical Gravitational Lensing Experiment, Udalski et al. 1992), Kepler (Borucki et al. 2003), and TESS (Transiting Exoplanets Survey Satellite, Ricker et al. 2015) have enormously contributed to the discovery of a large fraction of (candidate) binary systems in the nuclei of PNe. Numerous papers that have been published recently focus on this topic with encouraging results, significantly increasing the binary population (see, e.g. Miszalski et al. 2009; De Marco et al. 2015; Jacoby et al. 2021; Aller et al. 2020). However, the fraction of known binary central stars in the whole population of PNe remains very low, representing less than 1% of the whole sample of known Galactic PNe, while unbiased surveys indicate that the true fraction should be ∼20% or more (Jones & Boffin 2017; Jacoby et al. 2021). The need to continue this search for binary central stars is, therefore, essential in order to ultimately obtain a statistically significant and unbiased sample with which to definitively understand the role of binarity in the formation and evolution of PNe (as well as related phenomena).

In this paper, we continue the analysis carried out in Aller et al. (2020), in which we investigated the variability of the eight central stars of PNe observed in Cycle 1 of TESS for indications of binarity. Here, we extend this study to cover the first four years of TESS operations (Cycles 1 to 4) in order to find new short-period binary central star candidates in the whole sample of galactic PNe. In this search, as well as some already known binary central stars, we also identified new binary candidates that present modulations in their TESS light curves compatible with the presence of companion stars. These modulations are likely due to irradiation effects (i.e. the reflection of the light of the central star on the companion’s surface), ellipsoidal modulations (consequence of tidal forces with the companion), and Doppler beaming (relativistic effects due to the orbital motion of a binary system with respect to the observer). Previous works have used these effects to detect binary companions (see, for instance, Santander-García et al. 2015; De Marco et al. 2015; Maxted et al. 2002; Faigler & Mazeh 2011) and even planets (e.g. Shporer et al. 2011; Mazeh et al. 2012; Lillo-Box et al. 2014; Millholland & Laughlin 2017). All these physical effects open the door to the detection of close binaries without requiring the detection of eclipses or the need for radial velocity monitoring. Also, TESS represents a valuable opportunity to revisit close binary CSPNe already known by improving the determination of their properties, especially in those eclipsing systems.

The paper is organised as follows. In Section 2, we present our target sample and briefly introduce the TESS datasets. The light curve modelling process is explained in Section 3 and Section 4 is dedicated to the analysis of the results, discussing the detected variability in the light curves and the periodicities found, as well as distinguishing whether they are already known binaries, new binary candidates or stars with weak evidence or no variability at all. In Section 5, we discuss the origin of the variability in each particular case and, finally, we close the paper in Section 6 with some short final remarks and main conclusions.

2 Target selection and observations

2.1 The PN sample in TESS

TESS was launched in 2018 and was initially designed as a two- year mission, although it is now conducting its second-extension. During the first year of observations, TESS completed Cycle 1 in a total of 13 sectors and provided two-minute cadence light curves with a baseline of ∼27 days for eight central stars of PNe which were analysed in detail in Aller et al. (2020). In July 2019, TESS started its Cycle 2, observing the northern ecliptic hemisphere in another 13 additional sectors. Since then, every year TESS completes a new cycle. In this paper, we have analysed all the CSPNe that have been observed with two-minute cadence until September 2022, that is, Cycles 1, 2, 3, and 4. They are listed in Table A.1, along with their most relevant information as PNG designation, common name, TESS input catalogue (TIC) designation, equatorial coordinates, TESS magnitude, the sector or sectors of the observation and the PN status found in the HASH database¹ (Bojičić et al. 2017). In total, we analysed the available two-minute cadence light curves for 62 central stars of true (T), possible (P) or likely (L) PNe, including the eight previously presented in Aller et al. (2020) from Cycle 1. We note that almost half of the targets do not have the TIC designation linked to SIMBAD, so we have carefully inspected all the target pixel files (TPFs) and targets, to ensure that the light curves correspond to the target itself.

2.2 TESS photometric data

We retrieved the light curves provided by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016) from the Mikulski Archive for Space Telescopes MAST²). We used the PDCSAP (Pre-search Data Conditioning Simple Aperture Photometry) flux for the analysis, only removing those data points with a non-zero ’quality’ flag. For a more detailed description of the TESS data, we point the interested reader to Aller et al. (2020). When more than one sector was available in the MAST archive for a particular target, we combined all the individual light curves to produce a single and continuous multi-sector light curve. This is done with the stitch method in the Lightkurve package, which concatenates and normalises all the sectors.

In order to check for possible blends and contaminating sources in the apertures, we plotted the TPF of each central star in the sample using tpfplotter³ (Aller et al. 2020). This check of the TPFs is especially important in the most crowded fields, since the relatively large pixel size of the TESS CCDs (21 arcseconds per pixel) can lead to strong photometric contamination from nearby sources. Figure 1 shows an example TPF for three different degrees of contamination. In all cases, the central star is marked with a white cross and the red circles represent Gaia sources from the DR3 catalogue (Gaia Collaboration 2021), scaled by magnitude contrast against the target source. We have over-plotted all the sources with a magnitude contrast up to ∆m = 6 (i.e. six magnitudes fainter than our target, which corresponds to ∼ 0.4% of contamination if entirely inside the aperture). The aperture mask used by the pipeline to extract the photometry has also been plotted over the TPF.

As shown in the figure, we find very crowded fields with severe contamination (larger than 10%), fields with little contamination and cases with negligible contamination where no other star within a contrast magnitude of ∆m = 6 is identified by Gaia inside the TESS aperture. We note that this is only a first selection criterion that allows us to avoid critically contaminated cases where the attribution of the true source of the variability will not be possible with the current dataset and tools. For targets passing this criterion, we perform an analysis using the TESS_localize code to unveil the origin of the photometric variations and account for the possible contamination of stars outside the aperture. Figure B.1 shows the TPFs of the first sector observed for each target.

For each individual central star and sector, we have calculated the contribution of the field stars to the total flux in the aperture mask. Out of the 62 targets, we found severe contamination for 28 stars. We decided to exclude these targets from further analysis, since the probability that any potential variability can come from any of the nearby stars is too high. Only for specific contaminated objects, we carried out additional studies. In particular, Abell 63 (because it is a well-known eclipsing binary central star in an extremely crowded field), and also three objects that we already analysed in Aller et al. (2020) from Cycle 1: NGC 246, NGC 2867, and NGC 5189 (see Sect. 4). In addition, we found 16 cases with minimal (≤10% but not null) contamination and 18 central stars where there was no appreciable contamination from nearby sources in the aperture mask. The level of contamination (negligible, minimal or severe) is marked in Table A.1 with a, b or c, respectively. We subsequently discuss these cases in detail in Section 4.

Once this contamination cleaning step is done, each single multi-sector light curve was inspected by performing a Lomb– Scargle analysis in order to identify possible periodic variability in the data as a first look.

At this point, it is important to note that apart from the possible contamination of other stars in the field, we sometimes have bright PNe around our targets of interest, whose light might dilute the amplitude of the detected variability and hence, have an impact on the derived physical parameters. Quantifying this PN contribution is out of the scope of this study, which is to identify new binary candidates in the sample of central stars of PNe.

Fig. 1 Target pixel files (TPFs) of three central stars in the sample (marked with white crosses) obtained with TPFPLOTTER. The red circles are sources from the Gaia DR3 catalogue in the field with scaled magnitudes (see legend). The aperture mask used by the pipeline to extract the photometry is also marked. The pixel scale is 21 arcsec pixel−1 . Left, middle and right panels represent examples of fields with negligible, minimal or severe contamination, respectively (see Sect. 2.2 for more details).

3 Methodology: Data analysis and photometric modelling

3.1 Parametric approach

In a first approach, similar to the analysis performed in our previous work in Cycle 1 (Aller et al. 2020) and other works using the same detection technique (e.g. Millholland & Laughlin 2017; Lillo-Box et al. 2016), we model the TESS light curves by using simple sinusoidal functions for the three effects (ellipsoidal, Doppler beaming and reflection). This corresponds to $${\left(\frac{\text{Δ}F}{F}\right)}_{\text{ellip }}=-A_{\text{ellip }}\cos (2\theta ),$$(1) $${\left(\frac{\text{Δ}F}{F}\right)}_{\text{beam }}=A_{\text{beam }}\sin (\theta )\text{, and }$$(2) $${\left(\frac{\text{Δ}F}{F}\right)}_{\text{ref}}=-A_{\text{ref}}\cos (\theta ),$$(3)

with θ = 2π(t − T₀)/P, where T₀ is the time of inferior conjunction of the companion star, P is the orbital period, and A_ell , A_beam, and A_ref are the amplitudes of the ellipsoidal, Doppler beaming, and reflection effects, respectively. The simple model thus corresponds to: $$\begin{array}{l}{\left(\frac{\text{Δ}F}{F}\right)}_{3\text{E}}=Z_0+{\left(\frac{\text{Δ}F}{F}\right)}_{\text{ellip }}+{\left(\frac{\text{Δ}F}{F}\right)}_{\text{beam }}+{\left(\frac{\text{Δ}F}{F}\right)}_{\text{ref }}\\ =Z_0-A_{\text{ellip }}\cos (2\theta )+A_{\text{beam }}\sin (\theta )-A_{\text{ref }}\cos (\theta ).\end{array}$$(4)

As in most cases only one of the effects might be present, we also test a simpler model with only one sinusoidal function, namely: $${\left(\frac{\text{Δ}F}{F}\right)}_{1\text{E}}=Z_0+{\left(\frac{\text{Δ}F}{F}\right)}_{1\text{E}}=Z_0+A_{1\text{E}}\cos (\theta ).$$(5)

Based on these equations, we have seven parameters in the case of the 3-effects (3E, Eq. (4)), namely the photometric level (Z₀), the three amplitudes (A_ellip, A_ref, A_beam), the time of inferior conjunction (T₀ ), the orbital period P, and the photometric jitter (a term that we add in quadrature to the photometric uncertainties to account for other unknown systematics). For the 1-effect model (1E, Eq. (5)), we only have four parameters as only one amplitude is needed (A_1E).

We acknowledge the fact that this purely sinusoidal model is only an approximation to a physically motivated model. For instance, gravity darkening and geometrical effects produce non- sinusoidal ellipsoidal modulations. However, this approximation is sufficient for the exploratory goals of this work.

We explore the parameter space through a Monte Carlo Markov Chain sampler using the emcee algorithm (Foreman-Mackey et al. 2013). In particular, we set uninformative (uniform) priors for the orbital period between 0.1 and 10 days and for the time of inferior conjunction between the first date of observation and 10 days afterwards. In general terms, we allow the amplitudes to vary from 0 to twice the maximum peak- to-peak value of each light curve. The photometric level is also left with a uniform prior 𝒰(0.9,1.1), as well as the jitter term that we sample in the range 𝒰(0,100) in ppm.

In total, three models are tested for each of the targets studied in this paper: a flat model (FL) with only the photometric level and jitter as parameters, a 1-effect model (1E), and the 3-effects model (3E). For each of them, we sample the posterior distribution of the parameters using a number of walkers equal to four times the number of parameters of the model, and a total of 100 000 steps per walker. We check the convergence of the chains by ensuring that the length of the chain is at least 30 times the autocorrelation time (Foreman-Mackey et al. 2013). We estimate the Bayesian evidence (ln 𝒵_i) of each model using the perrakis⁴ code (Díaz et al. 2016). This metric is then used to select the simplest model that best represents the data. For a complex model to be selected over a simpler one (i.e. the 1E over the FL or the 3E over the 1E), we require the Bayesian evidence of the more complex model to be larger than the evidence of the simpler model by 6 in logarithmic space (i.e. 𝓑 = Δ In 𝒵 > 6, Trotta 2008), corresponding to strong evidence in favour of the more complex model. Otherwise, the simpler model is selected.

3.2 Identification of the variability source

As mentioned in Sect. 2.2, the large plate scale of TESS (21 arcsec pixel⁻¹ ) makes it particularly difficult to undoubtedly assign the source of the variability to one specific target (e.g. Lillo-Box et al. 2024). Essentially, the high frequency resolution of TESS is at the expense of its low spatial resolution. Even in those stars in which the aperture mask has little contamination from nearby sources, the analysis is delicate, since other stars outside the aperture can contributed to the measured flux. Recently, Higgins & Bell (2023) developed a method to localise the origin of variability on the sky to better than one fifth of a pixel given a measured frequency (or frequencies). Basically, the method can resolve the variable source in frequency space for each pixel. The authors showed that even stars more than three pixels outside the aperture, can produce significant contamination in the extracted light curves.

Once the MCMC analysis presented in Sect. 3.1 is performed and the periodicities in the TESS light curves are measured, we can apply this methodology. To this end we use the implementation presented in the open-source Python package TESS_localize⁵. We apply this algorithm to our sample (both the already known binaries and those with newly detected variability) in order to identify the origin (i.e. location on the sky) of the periodicities found in this work. This analysis is presented in Sec. 5.2.

4 Results in the TESS dataset

The results of the variability analysis described in Section 3 are summarised in Table 1. The column labelled as ‘Best model’ shows the preferred solution of the modelling process based on the Occam’s razor quantified by the comparison of the log of the Bayesian evidence between two models (𝓑_ij, where “i” and “j” represent each of the two models). This is, we compare the evidence of the 1E (1-effect) model against the flat (FL) model (𝓑10), the evidence of the 3E model against the flat model (B30) and the evidence of the 3E model against the 1E model (𝓑₃₁). Out of 38 targets, a total of 26 show modulations compatible with one or more effects. These central stars are described in Sections 4.1 and 4.2, depending on whether they are already known binaries or new binary candidates, respectively. For the remaining 12 central stars, the flat model is the most plausible, that is, there are no statistically significant modulations in the light curve down to the sensitivity of the data. These central stars are briefly discussed in Sect. 4.3.

4.1 Already reported binaries from the literature

Eight of the targets in the sample are either already known binary systems (confirmed through radial velocity observations) or candidates reported in the literature by using other techniques that still await confirmation. We do not include here our candidates from Cycle 1 of TESS presented in Aller et al. (2020), which will be discussed in Sect. 4.2 with the rest of the sample. The eight already known binaries are: Abell 63, Abell 46, AMU 1, DS 1, Abell 30, NGC 2392, NGC 5189, and LoTr 5. The first seven are short-period binaries and we are able to recover the expected variability in the TESS data. In all the cases except one (Abell 63), the contamination in the aperture mask from nearby stars is less than 10%. The remaining system is the long-period binary LoTr 5. In this case, we detect the variability associated with the rotation of one of the stars in the binary nucleus but obviously not the variability associated to the orbital period of the binary (we note that the central star of LoTr 5 has one of the longest orbital periods − P∼2700 days − of all the known central stars of PN, Jones et al. 2017). In the following, we revisit and briefly discuss each system, and present our new analyses. The results are summarised in Table 2.

Abell 63 and Abell 46. The eclipsing binary nuclei of PNe Abell 63 and Abell 46 have been extensively studied by several authors in the past years. Their nuclei, named UU Sge and V477 Lyr, respectively, consist of a primary hot white dwarf or subdwarf and a low-mass star. Their orbital periods are 0.465 days for UU Sge (Bell et al. 1994) and 0.472 days for V477 Lyr (Pollacco & Bell 1994).

Both systems have been the subject of detailed modelling in the literature (e.g. Afşar & Ibanoğlu 2008), in order to constrain the physical properties of the binary components (masses, temperatures and radii). In the case of Abell 46, for example, this has resulted in a potentially post-RGB scenario for the central star based on the low derived luminosity for the primary component. However, the derived temperature and surface gravity do not lie on post-RGB evolutionary tracks for the derived mass (instead being consistent with a much lower mass; Jones et al. 2022). Our photometric modelling does not take into account eclipses (see, for instance the TESS light curve from Abell 63 in Fig. 2) and, therefore, the analysis of these light curves cannot be performed by using the methodology explained in Sect. 3.

We, therefore, decided to model the system using a version of the PHOEBE2 code (Prša et al. 2016; Horvat et al. 2018; Jones et al. 2020; Conroy et al. 2020), adapted to incorporate model atmospheres calculated by Reindl et al. (2016, 2023) using the Tübingen Model Atmosphere Package (TMAP; Rauch & Deetjen 2003; Werner et al. 2003, 2012), and previously used to model the binary central star of the PN Ou 5 (Jones et al. 2022). This new approach also has the benefit of using more appropriate TMAP atmospheres for the primary star, where previous models either assumed Kurucz (1993, more suitable for main sequence stars) or black body atmospheres – both of which can lead to important differences in the derived temperature (potentially leading to the aforementioned mismatch with evolutionary tracks).

Preliminary explorations of the parameter space, fixing the stellar masses to those from previous modelling efforts (as they are based primarily on the observed radial velocities), indicated that a solution consistent with evolutionary tracks should be possible. As such, we ran a MCMC sampling of the binary parameters, allowing the primary temperature to vary but forcing its radius to lie on the appropriate evolutionary track of Miller Bertolami (2016), while leaving the temperature and radius of the companion free. The resulting best fit is shown in Fig. 3, and the best-fitting parameters are listed in Table 3. Ultimately, a more detailed fit is required (taking into account multi-band groundbased photometry as well as directly fitting the observed radial velocities – both of which are beyond the scope of this work) before claiming that the discrepancy between parameters from evolutionary tracks and binary modelling has been resolved or that the new parameters (and their associated uncertainties) are more reliable than those in the literature. Nevertheless, the fit presented here highlights that the high quality data from TESS combined with more advanced modelling tools has the potential to improve the derived parameters of post-CE binary central stars of PNe.

DS 1. In the nucleus of this PN there is a well-known double-lined spectroscopic binary (named KV Vel or LSS 2018) with an orbital period of ~0.357 days (Drilling 1985; Kilkenny et al. 1988). It is a post-common envelope system with an unusual and extremely strong reflection effect, with an amplitude of 0.55 mag in V (Hilditch et al. 1996). The TESS light curve clearly shows the same variability with a prominent signal at the same periodicity as the published orbital period. Preliminary fitting of the light curve using PHOEBE2 (following the scheme outlined above for the central star of Abell 46) indicates that the data are entirely consistent with the parameters derived by Hilditch et al. (1996), albeit with a slightly lower albedo for the secondary (0.55±0.01 c.f. 0.6).

Abell 30. Jacoby et al. (2020) reported the presence of light curve brightness variations in the K2 mission (Howell et al. 2014) data of the central star of this born again PN. Although it has not been confirmed through radial velocity observations, the authors concluded that these variations in the light curve were highly suggestive of a binary central system with a period of ~ 1.060 days. After analysing other possible physical processes, Jacoby et al. (2020) proposed the irradiation of a cooler companion as the most likely origin for the observed photometric variability in this system, and discarded other possible effects as Doppler beaming and ellipsoidal modulations. We find the same variability in the three sectors of the TESS light curves. Figure B.2 shows the Lomb-Scargle periodogram, where the False Alarm Probabilities (FAPs) at 10%, 1% and 0.1% are also indicated with grey horizontal dotted lines. The phase-folded light curve with the period derived in the fitting is plotted in the corresponding panel of Fig. 3. The one sinusoidal model appears as the preferred solution, showing a large amplitude of 14 ppt (see Table 2). However, we note that the evidence of a 3-effects model is almost the same as the simple model (𝐵₃₁ = −0.2), with an ellipsoidal amplitude of 0.95 ±0.25 ppt. Indeed, these ellipsoidal modulations can be easily recognised in the residuals of the one sinusoidal model in Fig. 3. However, although the time span covered by the three sectors of TESS (44, 45 and 46) is similar to that from the K2 mission, the precision of the TESS data is clearly worse than that from its predecessor and, therefore, the TESS light curves do not provide new information about the properties of the system.

AMU 1. The central star of this multipolar PN (Aller et al. 2013) was found to be a binary with a short period of ~2.928 days by De Marco et al. (2015). The Kepler light curve showed a photometric sinusoidal variability with a very low amplitude (0.73 ppt), consistent with relativistic beaming effects according to those authors. In the TESS light curves, we also detect a signal at the same periodicity, although the confidence level is much worse (approaching the False Alarm Probability at 10% in some sectors) because of the lower photometric precision and shorter baseline of TESS compared to Kepler. Though the results from the TESS light curve analysis show that the evidence for the flat model is larger, we can recover the 1-effect model with similar evidence discarding the sector 15, which has a large dispersion in the data. In this way, the MCMC chains converge to the known orbital period, although still with insufficient significance (𝓑₁₀ = −0.6, see Figure 4). In summary, the TESS data do not provide any new information beyond what Kepler already showed.

NGC 2392. The central star of this well-known PN was discovered to be a single-lined spectroscopic binary by Miszalski et al. (2019) after a radial velocity monitoring campaign.

According to the authors, the binary system appears to be a double degenerate system with an orbital period of ∼1.9 days. The periodogram of the TESS light curve shows a forest of significant periodicities above the 0.1% FAP (see Figure B.2). Among all these signals, we can identify two pairs of peaks more prominent at around ∼4.2 and ∼2.1 days, but none at the exact periodicity of 1.9 d found by Miszalski et al. (2019). We note, however, that the peak at 2.1 d is also present in the data by Miszalski et al. (2019), although the authors chose to analyse the 1.9 d because it was the highest in their RV periodogram. For that reason, we decided to re-analyse the radial velocity data presented in Miszalski et al. (2019). Our independent analysis using the same MCMC principles explained for the analysis of the TESS light curves but now applied to the radial velocity model, shows that the data converges to an orbital period of 2.10842273 ± 0.000085 days, which is actually significant in the RV periodogram (see upper panel of Fig. 5) and slightly longer than that derived by the authors with the same data. This result is in strong agreement with what we observe in the TESS light curve, and allows us to conclude that the 4.2 days periodicity seen in the light curve is, actually, an alias of the real period. This is a clear example of the importance of confirming the orbital period by more than one method. Taking this result into account, we forced the MCMC analysis to converge to the 2.1 days periodicity and obtained that the best solution is a 3E model.

NGC 5189. The central star of this complex, quadrupolar PN was already analysed in our previous work (Aller et al. 2020). Now, we re-analyse it taking into account new TESS data as well as the previous data from Cycle 1, drawing similar conclusions. After a radial velocity monitoring campaign, Manick et al. (2015) identified the periodicity found at 4.04 days as the orbital period of the binary. The periodogram of the TESS light curve (sectors 11 and 38) show several prominent peaks above the 0.1 % FAP (the most significant at 0.125, 0.858 and 1.71 d, see Figure B.2) but none of them around the 4.04 days period previously reported. The result of MCMC analysis is, therefore, very similar to those presented in Aller et al. (2020). A 3E model is the most plausible for the 1.71 d periodicity (the peak at 0.858 d is certainly an alias of this one) and the amplitudes are listed in Table 2. However, the resulted Doppler beaming amplitude is quite large and clearly not physically feasible for this type of object. This is a case with a lot of contamination from nearby stars in the photometric aperture so we cannot reach any firm conclusions from the TESS light curve. In fact, it is very likely that the origin of the periodicities we see in the TESS light curves does not correspond to our central star (see Section. 5.2 for more detail).

LoTr 5. This PN is known for having one of the longest period binary central star (P∼2717 days; Van Winckel et al. 2014; Jones et al. 2017) so far. The system consists of a hot star and a barium star with a rotation period of 5.95 days. A possible third component in the system was discussed by Aller et al. (2018) on the basis of new radial velocity observations, although without firm conclusions. The rotation period is clearly recovered in the TESS dataset (see Fig. 3). In the MCMC analysis, the reflection and the beaming effects are both degenerate, and the ellipsoidal variations go clearly to zero. Therefore, the 1E model is the only one possible, although there are clearly more effects in the light curve (as, for example, oscillation signals). These other effects certainly show up when subtracting the rotation signal from the light curve. We made this by using the detrending methods implemented in the WOTAN package (Hippke et al. 2019). Figure 6 shows the resulting periodogram of LoTr 5 after this detrending process. A more detailed analysis of these frequencies would be desirable to obtain precise information on the physical parameters of the pair, but this is out of the scope of this paper.

Fig. 2 Periodogram (upper panel), TESS photometric time series for sectors 40, 41, and 54 (middle panel), and phase-folded combined light curve (bottom panel) of Abell 63. Two different bin sizes are shown with black bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) circles.

Fig. 3 Phase-folded light curves of the already known binaries, with the best-fitting model overlaid. Results from the PHOEBE2 analyses are overplotted in pink and from the MCMC analyses in blue (1E solution) and red (3E solution). Grey symbols correspond to a bin size of 0.01 in phase (i.e. 100 datapoints) while the black symbols correspond to bin sizes of 0.05 in phase (i.e. 20 datapoints). The period derived from the fitting can be found in Table 2 and the corresponding periodograms are shown in Fig. B.2). We note that the period of NGC 5189 does not correspond with the orbital period published in the literature for this binary.

Fig. 4 Periodogram (upper panel), time-flux light curves for several sectors (middle panel), and phase-folded single, multi-sector light curve (bottom panel) of AMU 1. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting can be found in Table 2.

4.2 New binary candidates from TESS

In this section, we describe the results from the analysis of those central stars with negligible or minimal contamination in the TPFs and showing variability in the light curves compatible with the presence of (previously unknown) binary systems. In addition, we also revisit those binary candidates already analysed in Cycle 1 (Aller et al. 2020), even if they are found in severely contaminated fields. In total, we analyse in this section 18 central stars. The results from the analysis are summarised in Tables 4 and 5, depending on whether the best fit in the MCMC analysis is a 1E or a 3E solution, respectively.

We find that the variability of 12 central stars is better explained based on the current data with a simple sinusoidal model (see Table 4). Three of them (NGC 7293, NGC 2867 and RWT 152) are indeed binary candidates from Cycle 1 (Aller et al. 2020) that we have re-analysed by taking into account the new observations from Cycles 2, 3 and 4. For two of these targets (NGC 2867 and NGC 7293), the new results are compatible to those presented in Aller et al. (2020). In the case of NGC 2867, and after analysing sectors 9, 10, 35 and 37 together, we detect a significant periodic variability at ∼ 4.6 days with an amplitude of 0.3 ppt, similar to those from our previous paper and thus reinforcing the binary status of this central star. It is also important to highlight that in this case the 3E model has a larger evidence than the 1E model with a log-evidence of 𝓑₃₁ = +3.4 and the same periodicity. As shown in the corresponding panel of Fig. A.1, two bumps are appreciated at around phases 0.25 and 0.75, which is the reason for the 3E model having such a large evidence. However, the current dataset is still insufficient to prefer the more complex 3E model. It is also important to emphasise that the photometric aperture of NGC 2867 is quite crowded with other stars, so the origin of the variability is not entirely clear (see also Section 5.2). The results from the analysis of NGC 7293 (the well-known Helix Nebula) adding the new sectors 28 and 42, are also quite similar to those obtained in Aller et al. (2020). A clear photometric variability is identified with a periodicity of ∼2.8 days and amplitude of 1.5 ppt. Although, in this case, there is no contamination from nearby sources in the aperture, the origin of the variability may be also questionable (see Section 5.2). Only for the case of RWT 152, the new dataset from sector 34 leads to different conclusions than in Cycle 1. Now, apart from the forest of significant periodicities that we already saw in sector 7, a new clear signal at ∼6 days shows up in the periodogram (see corresponding panel in Fig. B.3), and the MCMC converges to that period by using a single sinusoidal model. As explained in Aller et al. (2020), some of the found periodic variabilities can be caused by the presence of spots. In any case, dedicated observations are required to confirm or discard the presence of a companion star.

Another nine central stars in the sample present variability compatible with a 1E model. They are: PG 1520+525, PNG 136.7+61.9, Pa 165, WPS 28, WPS 54, Fr 1-4, NGC 7094, NGC 2371 and Hen 3-1863. Fig. A.1 shows the phase-folded light curves with the periods derived in the fitting (see also Table 4). Objects with 1E model solution are plotted in blue. We briefly discuss here some key aspects of them as follows:

PG 1520+525. This is the only “classical” round (with a ring appearance, Jacoby & van de Steene 1995) PN in the sample showing variability in the TESS light curve. We find a periodicity of ∼5.2 d with an amplitude of ∼4.13 ppt, compatible with irradiation on a close companion. Although there are no other stars identified by Gaia inside the TESS aperture within ∆m = 6 magnitudes, the location of the variability is not clearly associated with the central star (see Sect. 5.2), so dedicated follow-up is desirable to confirm the origin of the detected variability.

PNG 136.7+61.9. This extremely faint PN⁶ with elliptical appearance, has a DAO (hydrogen and helium) white dwarf central star, suggested as a post-RGB candidate by Reindl et al. (2023). The authors did not find significant light curve variations in the TESS data. We find a 25 ppt variability amplitude in the only sector available in TESS that we could attribute to the irradiation of a cool companion.

Pa 165. Although the 1E solution is the most likely for this central star with an irregular PN, the 3E model has also large evidence (𝓑₃₁ = −3.6). The low amplitude of the variability (0.65 ppt) is compatible with either irradiation or relativistic effects.

WPS 28. The light curve of this white dwarf (with a possible PN around it according to the HASH PN database) exhibits a clear sinusoidal variability with a period of ∼2.8 days. The large temporal coverage of the TESS data, with a total of 25 sectors from Cycles 1 to 4, provides an exquisite precision. Interestingly, the phase-folded light curve with the period corresponding to this variability also shows a small dimming at the minimum of the relative flux (see corresponding panel in Fig. A.1). Although we note that this dimming is not yet statistically significant, this is the expected location (within the uncertainties), of an eclipse when the variability is due to the reflection effect and the orbital inclination is such that the companion passes between the stellar disk and our line of sight. In order to further test the source of the variability and this tentative eclipse signal, we use a modified version of the TESS positional probability software (tpp⁷) described in Hadjigeorghiou & Armstrong (2023) and Lillo-Box et al. (2024) to compute the probability of each Gaia-detected source in the field around the target as being the true source of the transiting signal. By doing so, we find that the star with the largest probability of being the transit host is actually WPS 28, although still with a 32% probability, which makes highly interesting the follow-up of this target. However, their TPFs show severe contamination from other stars in all the sectors except one (sector 21, which has minimal contamination). Additional observations are thus critical to ensure the source of the variability (but see also Sect. 5.2).

WPS 54. This star was classified as a RS CVn type star by Drake et al. (2014) based on their Catalina Sky Survey light curve. The RS CVn variable type consists of binary systems with high chromospheric activity FGK companions. After analysing the Catalina light curve, Werner et al. (2019) found a variability with a periodicity of 3.45 days, that the authors interpreted as the white dwarf’s rotation period. The TESS light curve actually shows this periodicity in the periodogram (see Fig. 7). Although we acknowledge that the variations in the light curve may be due to multiple spots on the star, we still try our modelling process. The 1E model on this periodicity produces large residuals, suggestive of a wrong model (see Fig. 7, upper panel). The 3E model on this periodicity still produce large residuals with the presence of a phase-shifted signal at half of the periodicity (1.7 d; see also Fig. 7, upper panel). We then decided to simultaneously model a one sinusoidal model for the 3.45 d signal (assuming it as the rotation period of the white dwarf) and 1E or 3E model to the 1.7 d periodicity. By doing so, we find that the residuals of both 1E and 3E model are largely reduced (see Fig. 7, bottom panel), with the 1E model being preferred against the 3E model. Still, this result may point to either irradiation with 1.7 d orbital period or perfectly symmetric ellipsoidal effect with a 3.45 d period. Both scenarios are thus very interesting because any of them suggest a binary companion with a coupled orbital period commensurable with the rotation period of the star, suggestive of a coupling between the orbital and rotation frequencies.

Fr 1–4. The central star of this possible (according to the HASH database) PN is the hot subdwarf JL 102. The single TESS light curve of this target (produced by combining sectors 27 and 39) shows a clear variability with ∼6.4 d period. The corresponding signal in the periodogram is very strong and our analysis reveals a 1E model with ∼5 ppt amplitude, which is compatible with irradiation on a stellar companion.

NGC 7094. This object was classified as a PN by Kohoutek (1963). NGC 7094 has a filamentary structure, with a shell deviating sphericity (Rauch 1999), and a PG1159-type central star (Löbling et al. 2019). The variability in the TESS light curve is very clear, with a strong peak at ∼ 4.3 d in the periodogram, which may be perfectly compatible with irradiation of a companion star. The contamination from other field stars inside the TESS aperture is negligible and NGC 7094 is one of the most promising candidates to reside in a binary system (see Sect. 5.2).

NGC 2371. This complex PN has been extensively studied in the literature (see, for example, Gómez-González et al. 2020; Ramos-Larios & Phillips 2012). It exhibits a high-excitation, multipolar morphology with a Wolf-Rayet-type as a central star. This central star is a well-known pulsator (Ciardullo & Bond 1996), whose TESS light curve was recently analysed by Córsico et al. (2021). The authors detected several pulsation periods for this star between 878.5 s and 1032.6 s (shorter than the pulsation frquencies previously reported by Ciardullo & Bond 1996, from the ground), but no evidence of binary signatures. On the contrary, we find several periodicities above the 0.1% FAP in the TESS periodogram, that populate a wide region at lower frequencies (most of them between 0.3 and 4 days, see Fig.B.3). Our analysis converges to the one at ~2.3 days, showing a low-amplitude (~0.8 ppt) sinusoidal variability.

Hen 3-1863. The light curve of this central star is highly dominated by a strong variability, which appears to have its origin in the central star (see Sect. 5.2). The periodogram shows a strong a broad peak at ~4.3 d (see Fig. B.3) and our MCMC analysis converges to that solution. The phase-folded light curve with that periodicity shows a structure at phase ε = 0.5 and large amplitude high-frequency variations that are not explained with our sinusoidal model (see corresponding panel in Fig. A.1), indicating that additional effects are present in the data. The nature of these effects are out of the scope of this paper. However, the large scale variability at 4.3 days is clear and well represented by our model.

On the other hand, we found 6 targets showing a 3E solution as the most likely one. The phase-folded light curves with the periods derived in the fitting of each of them are shown in Fig. A.1, with the 3E model solution overplotted in red.

Among these, two (NGC 246 and PG 1034+001) are binary candidates from Cycle 1 (Aller et al. 2020) that we re-analyse by adding the new data from Cycles 2-4. The results for PG 1034+001 are practically the same as in Aller et al. (2020), with a strong periodicity at ~1.86 days and evident ellipsoidal modulations. We invite the reader to revise our previous paper for more details. The case of NGC 246 is quite different. The central star of this elliptical PN is a known hierarchical triple stellar system (Adam & Mugrauer 2014). In Aller et al. (2020), we found a simple sinusoidal model with a period of ~6.8 days as the best solution. With the new data (sector30), the MCMC analysis converges to 3E solution at ~5 days. However, we note that both signals may have the same origin, since the peak is quite broad in the periodogram (see Fig. B.3), which is dominated by the low frequencies (in the range of 500–700 µHz) associated with pulsations (Ciardullo & Bond 1996). Also, it is important to remember that this is one of the targets with substantial photometric contamination by other nearby stars in the TPF (see corresponding panel in Fig.B.1), which makes it extremely complicated to draw firm conclusions on the possible binary scenario.

Other 4 objects in the sample present variability compatible with a 3E modulation. They are briefly discussed as follows:

NGC 1501. The central star of this PN is a pulsating hydrogen-deficient, pre-white dwarf. Its TESS light curve was already analysed by Córsico et al. (2021) in order to search for pulsation frequencies and binary signatures. Although evidences for the latter scenario were not found, the authors presented asteroseismological analysis for the central star. We find several significant peaks in the periodogram, apart from the low-frequency signal analysed by Córsico et al. (2021), which may be compatible with the presence of another star in the system. Our analysis converges to the one at ~3.3 days by using a 3E solution, showing a low (~0.5 ppt) beaming amplitude that is compatible with relativistic effects in the system.

K1–16. As NGC 1501, this is also a pulsating hydrogen-deficient, pre-white dwarfs whose TESS light curves were already analysed by Córsico et al. (2021). Binary signatures were not found and their asteroseismological analysis was quite limited in this case because of the dramatic changes in the pulsations of this star. The periodogram shows a forest of significant periodicities (between 1 and 11 days, approximately, see corresponding panel in Fig. B.3). The MCMC analysis found the second most prominent peak (corresponding to ~ 3.9 days) as the most probable, with low amplitudes for the three effects (irradiation, ellipsoidal modulations, and Doppler beaming), although the beaming amplitude still is unrealistically large. Therefore, the origin of this variability is hard to interpret due to the large amount of different signals. Radial velocity observations are crucial to confirm the nature of this variability.

IC 2149. The central star of this elongated PN shows a high degree of variability in the TESS light curve, resulting in a large number of periodicities above the 0.1% FAP in the periodogram (see Fig. B.3). Given the large amplitude and high-frequency of these variations and the fact that we are interested in longer periodicities, we included in the modelling an extra jitter noise term to account for these short-period signals in a proper way. After doing so, the preferred solution for the MCMC fit is a 3E model with a period of 1.31 days. The derived Doppler beaming amplitude (A_beam = 1.01 ± 0.25 ppt) is indeed too large to be realistic and might still be driven partly by the high-frequency variability. However, it is important to note that its uncertainty is still large and might be compatible with 0 at 4σ.

IC 4593. As the case of the previous two PNe, the peri-odogram of this central star is dominated by several peaks above the 0.1% FAP. However, there are three of them (between 2 and 4 days) that predominate above the rest. The analysis converges to the one at 2.6 days, with a large reflection amplitude (A_ref ~12.8 ppt), a moderate A_beam of ~2.7 ppt and no Doppler beaming.

Fig. 5 Lomb-scargle periodogram (upper panel) and radial velocity analysis (bottom panel) of NGC 5189, based on the data presented in Manick et al. (2015).

Fig. 6 Periodogram of LoTr 5 after subtracting the rotation period with the WOTAN package (in blue and orange, with two different window sizes in the flatten module) and the original periodogram (in black).

Fig. 7 Modelling of the WPS 54 TESS light curve. Top panels: phase- folded light curve (top-left) with the periodicity corresponding to the largest peak in the periodogram (top-right panel), including the 1E (blue) and 3E (red) models and their corresponding residuals. The model is not satisfactory due to the presence of additional signals. Bottom panels: phase-folded light curve (bottom-left) with the 1.73 d periodicity after removing the sinusiondal patter corresponding to the 3.45 d (potentially the rotation period). The 1E (blue) and 3E (red) models are shown, as well as their corresponding residuals.

4.3 Targets with no variability or with weak evidence

We find 12 central stars (leaving aside AMU 1, already discussed in Section 4.1) in which the best solution is the flat model (see Table 1). In general terms, we can discard sinusoidal variabilities down to the TESS sensitivity for these targets. Although in some few cases like AMU 1 and Fr 2-21 the MCMC finds compatible solutions with sinusoidal patterns, the Bayesian evidence still favours the simpler flat model. For all these non-detections, we can use the posteriors of the MCMC analysis to provide an upper limit to the amplitude of potential variability of a hypothetical companion. We use the marginalised posterior distribution of the 1E model (A_1E) and use its 95% percentile as the upper limit. This is shown in Fig. B.5, where we display the corner plots corresponding to the dependency between the period and A_1E as well as the marginalised distributions and the 95% upper limit.

At this point, there are some cases that deserve special attention. The case of Abell 7 is intriguing since, in Cycle 1 (sector 5), we identified two clear significant (and independent) peaks in the periodogram at ∼2.6 and ∼3 days (Aller et al. 2020) with unclear origins. Interestingly, these peaks no longer appear in the new sector 32 data (see Fig. B.4), so that the possible binary scenario remains questionable. As a possibility, the detected variability in sector 5 may be due to other non-binary scenarios as, for example, starspots.

We would also like to mention the case of StDr 56⁸. It is a faint and filamentary nebula, in which central region Drech- sler and Strottner identified two possible white dwarfs, either of which could be the central star of the nebula. It is worthy to note that, although the flat model is the best solution in the MCMC fitting, the periodogram of the TESS light curve shows a very clear significant peak at ∼7 days (see Fig. 8, top panels). However, the uncertainty of the data is still large compared to the amplitude of the candidate signal and so the flat model is preferred. This is a promising candidate for subsequent follow-up.

Another good binary candidate is the central star of Fr 221. Although, as in the case of StDr 56, the best solution for this target is the flat model (see Table 1), the 1E and 3E models have comparably large evidences compared to the flat model (𝓑₁₀ = 4.5, 𝓑₃₀ = 5.0). The 1E solution is shown in Fig. 9, where we present the TESS time-flux light curve, the periodogram and the phase-folded light curve at the strongest peak (∼2.5 days). For this analysis, several photometric points present just before and after the TESS data downlink have been removed, since they are typically the result of instrumental effects. It is also worth mentioning that a possible eclipse located at the expected minimum of the phase folded light curve is also hinted. However, the current data do not allow its confirmation.

Finally, the central stars of PNe JnEr 1 and Abell 28 were found to have infrared excess compatible with a possible companion spectral type later than M4 by Douchin et al. (2015). The TESS light curves of these objects do not show any variability suggestive of a close binary scenario. Only in the case of Abell 28 do we find a significant signal at 0.11 days above the 0.5% FAP in the periodogram (see Fig. B.4), but the MCMC analysis does not converge to that peak. Similarly, De Marco et al. (2013) detected infrared excess in Abell 31 suggesting an M4V companion, in agreement with the works of Frew (2008) and Ciardullo et al. (1999). The high-precision photometry of the four sectors (34, 44, 45, 46) of TESS show no evidence of photometric variability at least in the time span covered by the light curves (see Fig. B.4). This is in accordance with the wide companion proposed by Ciardullo et al. (1999).

The rest of the targets which have no significant photometric variability in their light curves are HDW 7, Lo 1, EC 13290-1933, Sh 2-216, IC 5148/50 and Fr 2-46. Their periodograms, which do not show any peak above the FAPs, are also presented in Fig. B.4. Only for the case of Sh 2-216, the most prominent peak (at ∼0.2 days) could be a significant signal with more data.

Fig. 8 Periodogram (upper panel), time-flux light curve (middle panel), and phase-folded light curve (bottom panel) of StDr 56. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting can be found in Table 6.

Fig. 9 Periodogram (upper panel), time-flux light curve (middle panel), and phase-folded light curve (bottom panel) of Fr 2-21. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting can be found in Table 6.

5 Discussion of the variability

5.1 Physical reliability of the Doppler beaming amplitudes

It is important to note that, as the modelling process does not put meaningful constraints on the physical parameters, the solution from the MCMC analysis has to be checked for unfeasible solutions. Particularly, we discuss the plausibility of the relativistic effects (Doppler beaming) as responsible for the sinusoidal variability in both 1E and 3E solutions. This effect produces low- amplitude sinusoidal variations (Zucker et al. 2007) in the light curve, so in most of the cases we can directly reject this scenario by just checking for large amplitudes inferred from the MCMC analysis.

In the same fashion as in Aller et al. (2020), following the equations (2) and (3) from van Kerkwijk et al. (2010), we can estimate the radial velocity amplitude derived from the Doppler beaming amplitude obtained in the analysis (A_beam in the 3E solution and A_1E in the simple 1E solution under the assumption that the sinusoidal variability is due only to the Doppler beaming effect), and, therefore, estimate the binary mass function ( f ), defined as (assuming circular orbit): $$f=\frac{P_{\text{orb }}{v}_r^3}{2\pi G}=\frac{M_2^3{\sin }^{3}i}{{\left(M_1+M_2\right)}^2},$$(6)

where M₁ and M₂ are the masses of the two components, i is the inclination of the system, P_orb the orbital period and v_r the corresponding radial velocity semi-amplitude calculated from the measured Doppler beaming amplitude in the TESS light curve (Eq. (2) in Aller et al. 2020). By definition, the binary mass function will be always lower than M₂.

We derive this mass function for all the targets with 1E or 3E solutions after the analysis of their light curves. To simplify, and as there is no available information on the effective temperatures of all the central stars in the sample, we have used for the calculation two reasonable values for this parameter that approximately cover the whole possible range of it in this type of stars, that is T_eff = 40 000K, for the lower limit, and 150 000K for the upper one).

Figure 10 shows the calculated mass functions for 0.5, 1, 5 and 10 M_⊙ in the Amplitude vs P_orb plane. For visualisation purposes, cases with A_beam (in golden) or A_1E (in purple) greater than 6 are not shown in the plot (IC 4593, with A_beam ∼ 0.005 ppt, is also not included). As expected, large Doppler beaming amplitudes imply unrealistic masses for the hypothetical secondary component (assuming that the primary has a typical mass of a white dwarf or similar). Only for very short orbital periods (a few hours) would the solution be reasonable but we do not have orbital periods in that region. Thus, only a handful of central stars have amplitudes (A_beam or A_1E) compatible with relativistic effects. Subsequent follow-up efforts must take these unphysical scenarios into consideration.

5.2 Identification of the variability source

As explained in Sect. 2.2, we applied the TESS_localize method to our sample in order to identify the origin of the detected variability. The results are summarised in Table 7 where we list the name of the target in our sample, its TIC designation, whether or not the variability is attributed to the corresponding target (according to TESS_localize) and the most likely Gaia source which would be responsible for the variability, together with its probability.

In some cases, the result from the TESS_localize is not reliable according to its standards (see Higgins & Bell 2023 for further information), and we classify such targets as inconclusive. In these cases, the fit to the location of the source of the variability is not significant based on the height parameter of the code. According to the documentation of the algorithm, a proper fit is achieved if the relative uncertainty in the height is smaller than 20%, which corresponds to 5σ significance. In brief, only cases where the signal with the indicated frequency can be properly fit, can provide a good location for the variability. Otherwise, the solution is unreliable. We applied the algorithm to each individual sector of each target⁹, by using the optimal number of PCA components that TESS_localize calculates. For simplicity, Table 7 shows a mean of the likelihood calculated only with those sectors in which the height is ≤20%. If none of the sectors has a height ≤20%, the solution is unreliable.

The already known binary central stars in our sample are the best targets to test the method. The results from TESS_localize for these stars reveal that, as expected, the source of the periodicity found in their extracted TESS light curves match the corresponding positions on the sky of each specific central star. However, for NGC 5189, this is not the case. According to the TESS_localize results, the two periodicities detected in the TESS light curve (none of them coinciding with the published orbital period, detected through radial velocities, see above) have origins different from the central star. The photometric variability in TESS associated with the RV frequency (measured by Manick et al. 2015) is however too shallow to apply the TESS_localize test. Also, in the case of AMU 1, the signal of the variability is too weak to provide a reliable result, so TESS_localize does not recover the position of the central star as the origin of such variability.

The case of NGC 2392 is also of interest, since the variability at period 2.1 d (the one detected in the TESS light curve) matches with a 100% probability the position of the central star. Conversely, when asking the algorithm about the origin of the 1.9 days period previously reported by Miszalski et al. (2019), the probability that the origin of such periodicity being the central star is negligible. This result reinforces the 2.1 days periodicity as the true orbital period of this binary central star.

In addition, we find that for seven of our new binary candidates, the origin of the variability appears to coincide with the position of the central star (see Table 7). Thus, these targets have a very high probability to reside in binary systems and radial velocity follow-up would be desirable to confirm their binary nature and derive orbital and physical parameters. For one central star (NGC 7293), the detected variability is significantly attributed to another star in the field (Gaia DR3 6628875408232926336). Although it is a target with negligible contamination in the aperture mask, the probability that the measured variability corresponds to the mentioned field star is basically 100%, according to TESS_localize. This star is outside the aperture (more than 5 pixels far away from the central star) which highlights the importance of doing this kind of analysis. Thus, we can discard NGC 7293 as a potential binary candidate.

We would also like to mention the case of WPS 28. Out of 25 available sectors for this target, only three provide a reliable fit (with the relative uncertainty of the “height” parameter in these cases close to the 20%), so we think this result should be taken with caution.

For the remaining targets (8), the result of TESS_localize is unreliable because of the low significance of the detection. Thus, we cannot unequivocally confirm or discard the binary scenario for these targets and new dedicated observations are required.

Fig. 10 Amplitudes versus orbital periods measured in the TESS light curves. Doppler beaming amplitudes (from the 3E solution) are plotted in gold colour, while A1E (that could be due to irradiation either Doppler beaming), are in purple. Four different mass functions (0.5, 1, 5 and 10 M⊙) for two different effective temperatures (Teff = 40 and 150 kK) have been overplotted to see the position of the stars according to their corresponding effective temperature.

6 Final remarks and conclusions

As a continuation of the work presented in Aller et al. (2020), we analyse here the 2-minute cadence TESS light curves available in the MAST archive from the first four years of data from this mission (Cycles 1 to 4) of a sample of planetary nebulae nuclei in our Galaxy. In total, we analyse 62 central stars of (true, possible or likely) PNe in order to identify new binary candidates.

Figure 11 shows a summary of our results. Among the 62 central stars, 24 were set aside since they were in crowded fields, with large photometric contamination from nearby stars and, therefore, the origin of a hypothetically detected variability in such cases would be difficult to assess. The large size of the TESS pixel scale (21 arcsec pixel⁻¹ ) makes this first exploration essential. We then performed MCMC analyses to the remaining central star light curves (i.e. those with negligible or minimal contamination) in order to study the binary scenario, accounting for different physical effects: irradiation, ellipsoidal modulations and Doppler beaming.

Among the modelled sample, eight central stars are already known binary systems; another 18 present modulations in the light curves that may be compatible with binary signatures; and 12 central stars show weak or no significant variability.

For the already known binary systems, we analysed the new high-precision data from TESS, also using the radial velocity observations available in the literature for a handful of them.

This analysis allows us to provide new information in some of the systems, for example, refining the orbital parameters in the case of NGC 2392.

Additionally, we propose as binary candidates the 18 central stars showing sinusoidal variability in their light curves (four of them already analysed in Cycle 1). For all of them, we derived short orbital periods ranging from ∼1.3 to 6.5 days. With the goal of confirming the source of the variability, we used the code TESS_localize. This algorithm has been designed to identify, with an estimated likelihood, the location of the origin of a given frequency in the TESS target pixel file. Thus, after applying this method, we find that 7 central stars seem to be the origin of the detected periodicities, being our best candidates, while in three central stars the cause of the variability is likely another star in the field. For the remaining targets (8), TESS_localize provides unreliable fits so the results are inconclusive.

It is worth highlighting that the periods of the detected candidates do not match with the current period distribution of the known close binary central stars (see, for example, Jacoby et al. 2021). The candidates detected in our work have larger periods than the median of the distribution of published binaries. However, it is also true that most of them have been discovered through ground-based observations, that are biased towards short periods and large amplitudes. With the recent space surveys like Kepler, K2 and TESS, we are populating a region with lower amplitudes and longer periods, changing the period distribution we know so far. In any case, it is remarkable that in our sample, we do not find binary candidates with periods below ∼1 day. A priori, there is not any selection effect in the sample which may affect this finding, and it is true that other effects like stellar spots may be the answer to some of the variability encountered. To resolve this, more detailed analysis has to be done in each case to confirm the potential binarity and to constrain the parameters of the systems.

We note that among the binary systems (both already known and new candidates) only two PNe are classified as round in the HASH PN database, which represents less than 10%. They are Abell 30 and PG 1520+525. The first one, certainly has a spherical shell but with other structures like dots and cometary tails in the inner regions. On the other hand, PG 1520+525 is a more classical and perfectly round PNe. The rest, are PNe with morphologies that deviate from sphericity, such as bipolar, elliptical, and/or with other complex structures, that can be in many cases perfectly explained by the presence of a binary system in their nuclei.

On the contrary, if we attend to the morphology of the sample of the 12 central stars that do not show variability (i.e, they converge to a flat model in our MCMC analysis), the results are quite different. About ∼33% have round PNe and another ∼33% have non-defined morphological type according the HASH PN database. The remaining ∼33% have non-spherical shapes, such as elliptical or bipolar. The non-detected variability in such cases is by no means a confirmation that they are indeed single stars but instead, we provide a detection limit for a possible companion.

Definitively, TESS is providing a huge amount of high- precision photometric data that are enormously contributing to the study of a broad range of variable stars, including binary central stars of planetary nebula. The discovery and characterisation of more of these systems is crucial to further our limited understanding on the formation of the beautiful, varied and complex sample of PNe. The PLATO mission (Rauer et al. 2014) in the near future will also help in this way by providing long-term high-precision photometry for a large region of the sky.

Fig. 11 Results after the MCMC variability analysis of the total sample of central stars of planetary nebulae. Green dots indicate those targets with a good location of the variability in the central source after applying TESS_localize and red crosses indicate those targets with high probability that the variability is due to other star in the field. Red X marks indicate those objects with weak evidence which may become binary candidates with more data.

Data availability

Appendix B is available at https://zenodo.org/records/13370453

Acknowledgements

We thank Nicole Reindl and Mikkel Nørup Lund for useful discussions on the pulsations of LoTr 5, and Nicole also for her help in the selection of the targets and providing feedback on the manuscript. AA acknowledges support from Government of Comunidad Autónoma de Madrid (Spain) through postdoctoral grant ‘Atracción de Talento Investigador’ 2018- T2/TIC-11697. J.L.-B. was partly funded by grants LCF/BQ/PI20/11760023, Ramón y Cajal fellowship with code RYC2021-031640-I, by the Spanish MICIU/AEI/10.13039/501100011033 grant PID2019-107061GB-C61, and by the MICIU/AEI/10.13039/501100011033 and NextGenerationEU/PRTR grant CNS2023-144309. DJ acknowledges support from the Agencia Estatal de Investi- gación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) and the European Regional Development Fund (ERDF) with reference PID-2022-136653NA-I00 (DOI:10.13039/501100011033). DJ also acknowledges support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Inno- vación y Universidades (MCIU/AEI) and the European Union NextGenera- tionEU/PRTR with reference CNS2023-143910 (DOI:10.13039/501100011033). This research is part of the I+D+i project PID2019-105203GB-C21 founded by Spanish AEI (MICIU) grant 10.13039/501100011033. This research has made use of the SIMBAD database, operated at the CDS, Strasbourg (France), Aladin, NASA’s Astrophysics Data System Bibliographic Services, and the Spanish Virtual Observatory (http://svo.cab.inta-csic.es) supported by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 through grant PID2020-112949GB-I00. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. We acknowledge the use of pipelines at the TESS Science Office and at the TESS Science Processing Operations Center.

Appendix A Additional material

Fig. A.1 Results from the MCMC analysis for the proposed binary candidates. 1E solution is plotted in blue and 3E solution in red. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting in each case can be found in Tables 4 and 5.

References

All Tables

All Figures

Fig. 1 Target pixel files (TPFs) of three central stars in the sample (marked with white crosses) obtained with TPFPLOTTER. The red circles are sources from the Gaia DR3 catalogue in the field with scaled magnitudes (see legend). The aperture mask used by the pipeline to extract the photometry is also marked. The pixel scale is 21 arcsec pixel−1 . Left, middle and right panels represent examples of fields with negligible, minimal or severe contamination, respectively (see Sect. 2.2 for more details).

In the text

	Fig. 2 Periodogram (upper panel), TESS photometric time series for sectors 40, 41, and 54 (middle panel), and phase-folded combined light curve (bottom panel) of Abell 63. Two different bin sizes are shown with black bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) circles.
In the text

Fig. 3 Phase-folded light curves of the already known binaries, with the best-fitting model overlaid. Results from the PHOEBE2 analyses are overplotted in pink and from the MCMC analyses in blue (1E solution) and red (3E solution). Grey symbols correspond to a bin size of 0.01 in phase (i.e. 100 datapoints) while the black symbols correspond to bin sizes of 0.05 in phase (i.e. 20 datapoints). The period derived from the fitting can be found in Table 2 and the corresponding periodograms are shown in Fig. B.2). We note that the period of NGC 5189 does not correspond with the orbital period published in the literature for this binary.

In the text

	Fig. 4 Periodogram (upper panel), time-flux light curves for several sectors (middle panel), and phase-folded single, multi-sector light curve (bottom panel) of AMU 1. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting can be found in Table 2.
In the text

	Fig. 5 Lomb-scargle periodogram (upper panel) and radial velocity analysis (bottom panel) of NGC 5189, based on the data presented in Manick et al. (2015).
In the text

	Fig. 6 Periodogram of LoTr 5 after subtracting the rotation period with the WOTAN package (in blue and orange, with two different window sizes in the flatten module) and the original periodogram (in black).
In the text

Fig. 7 Modelling of the WPS 54 TESS light curve. Top panels: phase- folded light curve (top-left) with the periodicity corresponding to the largest peak in the periodogram (top-right panel), including the 1E (blue) and 3E (red) models and their corresponding residuals. The model is not satisfactory due to the presence of additional signals. Bottom panels: phase-folded light curve (bottom-left) with the 1.73 d periodicity after removing the sinusiondal patter corresponding to the 3.45 d (potentially the rotation period). The 1E (blue) and 3E (red) models are shown, as well as their corresponding residuals.

In the text

	Fig. 8 Periodogram (upper panel), time-flux light curve (middle panel), and phase-folded light curve (bottom panel) of StDr 56. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting can be found in Table 6.
In the text

	Fig. 9 Periodogram (upper panel), time-flux light curve (middle panel), and phase-folded light curve (bottom panel) of Fr 2-21. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting can be found in Table 6.
In the text

Fig. 10 Amplitudes versus orbital periods measured in the TESS light curves. Doppler beaming amplitudes (from the 3E solution) are plotted in gold colour, while A1E (that could be due to irradiation either Doppler beaming), are in purple. Four different mass functions (0.5, 1, 5 and 10 M⊙) for two different effective temperatures (Teff = 40 and 150 kK) have been overplotted to see the position of the stars according to their corresponding effective temperature.

In the text

Fig. 11 Results after the MCMC variability analysis of the total sample of central stars of planetary nebulae. Green dots indicate those targets with a good location of the variability in the central source after applying TESS_localize and red crosses indicate those targets with high probability that the variability is due to other star in the field. Red X marks indicate those objects with weak evidence which may become binary candidates with more data.

In the text

	Fig. A.1 Results from the MCMC analysis for the proposed binary candidates. 1E solution is plotted in blue and 3E solution in red. Two different bin sizes are shown with black (bin size of 0.01 in phase) and grey (bin size of 0.05 in phase) symbols. The period derived from the fitting in each case can be found in Tables 4 and 5.
In the text