© G. Kovacs et al. 2022

1 Introduction

The year 2005 marks the first direct detection of the light radiated by an extrasolar planet (Charbonneau et al. 2005; Deming et al. 2005). The observations were made by the Spitzer space telescope (Werner et al. 2004) at 4.5, 8, and 24 µm, which are not easily accessible by ground-based instruments. Although it was expected that a similar measurement in the near-infrared (NIR) could also be possible by ground-based 4m-class telescopes, two years passed until the first tentative observation of that kind was published (Snellen & Covino 2007). Since then, secondary eclipse (occultation) observations at the 2.2 µm (2MASS K – or Ks) band still remain in the realm of ground-based instruments, because of the lack of space instruments at this wavelength (e.g., Croll 2015; Zhou et al. 2015; Martioli et al. 2019). The 2MASS bands are especially suitable for the observation of hot extrasolar planets, because of the expected peak of the black-body flux in ~ 1–2 µm for temperatures between 1500 and 2000 K, that is, for the characteristic equilibrium temperatures of extrasolar planets (e.g., Alonso 2018).

Here we revisit WASP-5, an “ordinary” extrasolar planetary system, discovered by the SuperWASP collaboration (Anderson et al. 2008). The system harbors a single planet with a main sequence host akin to our Sun. So far, no other planets have been reported in the system, although there are contradictory results concerning the origin of the transit-time variation of planet b (i.e., Fukui et al. 2011; Hoyer et al. 2012). Based on the followup work of Gillon et al. (2009), the main system parameters are as follows: R_s/R_⊙ = 1.029, M_s/M_⊙ = 0.960, T_eff = 5700 K, a = 0.0267 AU, R_p/R_J = 1.087, and M_p/M_J = 1.58. These parameters imply an equilibrium temperature (assuming zero albedo and full heat redistribution) of 1740 K (Chen et al. 2014). The orbital period is 1.6284300 d, which was derived as part of the present study from the combination of the earlier epochs and those resulting from the analysis of the data from the TESS satellite (Ricker et al. 2015).

Occultation observations in the Ks band were previously carried out by Chen et al. (2014) and Zhou et al. (2015). Here, we combine these data with our unpublished observations made using the 6.5-m Walter Baade Telescope at the Las Campanas Observatory.

Our main goal is to increase the precision of the estimation of the occultation depth, which is an important ingredient for a more reliable model fitting, as most of the observations (including ours) sample the planetary spectra in only a few isolated bands. To constrain the atmosphere further, we use the recent data collected by the TESS satellite. We search for reflected light variation and occultation events. As a by-product of our analysis, we were also able to search for additional planets (but find none).

2 Data Sets

Two occultation light curves in the NIR Ks-band have been published so far on WASP-5 b. Chen et al. (2014) used the MPG/ESO 2.2 m telescope to observe the target in all three 2MASS bands. In spite of the substantial instrumental systematic errors, these authors clearly detected the event after applying corrections to counter positional and image quality dependencies. Zhou et al. (2015) performed a survey of seven hot Jupiters using the Anglo-Australian Telescope (AAT). Their survey also included WASP-5, yielding a long-duration coverage, allowing sufficient baseline in the eclipse modeling. Details of the observational settings and the methods used are described in the corresponding papers.

The third dataset comes from our single-night observations on 9 November, 2011 (UT). The four-chip camera of the FourStar infrared imager attached to the 6.5-m Walter Baade telescope was used to gather high cadence Ks images on the 10.9′ × 10.9′ field hosting WASP-5. An integration time of 4.4 s was used, yielding a ~7 s overall sampling interval. For better photometric accuracy, the telescope was slightly defocused, resulting in stellar images of ~10″ diameter. All images were taken in a simple staring mode, without dithering. Unfortunately, the sky was not photometric throughout the night because of intermittent clouds. This led to losing some 300 data points primarily after the ingress, affecting ~25% of the full observing run.

To obtain the photometric fluxes, we employed both the classical iraf¹ routines and those of the fitsh² package by Pál (2012). The two methods have led to very similar results, and so we decided to use our earlier reduction made by iraf.

First we performed the standard reduction steps of bias, dark, and flat corrections, including a treatment for the overall infrared sky emissivity variation by a nonlinear iterative multistep method, the nebulosity filtering algorithm³ of Irwin (2010). We then tested several aperture sizes to select the one that yielded the least scatter in the corresponding ensemble light curves (LCs). It turned out that nearly all apertures yield LCs of almost equal quality, with a slight preference towards mid-sized apertures. Finally, we selected the one with an aperture radius of 30 pixels (4.8″), and an outer annulus starting at a pixel radius of 40 and ending at 50 to assess the temporal background level.

In deriving the final ensemble LC (i.e., the target flux divided by the simple sum of the fluxes of the comparison stars), we decided not to use any comparison star from chips others than chip 1, which hosts the target. On this chip (see Fig. 1), we have two bright comparison stars (Nos. 2 and 3) and a fainter one (No. 4). We find that adding the fainter star slightly increases the noise⁴, and therefore we settled with the ensemble of the two brightest stars only.

Fig. 1 Field of the first chip of the FourStar infrared mosaic imager of the Baade telescope. North is to the left, east is to the bottom. Chips 2, 3, and 4 are located clockwise starting at the bottom of chip 1. The image size is 5.5′ × 5.5′. We used comparison stars Nos. 2 and 3 only.

3 Merging the Three Ks Light Curves

Before some of the peculiarities of the merging process are detailed, we describe the steps leading to the ensemble LC of the FourStar/Baade data (set-2 in Table 1).

3.1 The FourStar/Baade Light Curve

As mentioned, the most serious problem with the data is the temporal cloudiness during some part of the first half of the observation. The top three panels of Fig. 2 show the flux variation for the entire run, including the target and the two comparison stars. It is worth noting that for better visibility of the part of the flux variation that is dominated by the nonoutlying points, we limited the plots at the 4% flux drop. Several data points reach as much as 60–80% drops.

Although the comparison stars serve as an excellent diagnosis of the environmental origin of the harsh variations seen in the target, and the ensemble flux ratio cures most of the variation originating from the Earth’s atmosphere, we see in the bottom panel that the large drops in the flux could not be filtered out at the level required by the small signal we are searching for. Nevertheless, the binned LC strongly suggests the presence of an underlying occultation signal. We note that in constructing the binned LC, we used overlapping bin sets with a shift of half of the bin width. In this way, we can test the dependence of the binned LC on the bin distribution, which is an important piece of information when considering the sensitivity of any conclusion to be drawn from the binned LC, even if the conclusion is only preliminary.

In further processing the set-2 data, we observe the following: (1) There are outlier data points that are concentrated in a sufficiently broad section of the full time series, and therefore might seriously bias the derived eclipse parameters. (2) Likely because of differential extinction, albeit small, there is a significant downward trend in the ensemble LC. This should also be filtered out. (3) Finally, closer inspection of the ensemble LC in the bottom panel of Fig. 2 reveals that roughly in the middle of the cloudy period the flux suddenly jumped by a small fraction, enough to have a visible effect on the expected shallow eclipse. The most likely cause of this jump is the short-time change in the telescope pointing, leading to a sudden variation in the ensemble of pixels used in the flux evaluation. Leaving this jump in the ensemble LC would bias the occultation depth.

By following the principle of ‘least data massaging’, we proceeded as follows. For the treatment of outliers (issue 1), we robustly fitted⁵ eleventh-order polynomials to the fluxes of the target and the comparison stars. After the fit, we employed a 3σ-clipping for the outliers and replaced these items by the corresponding polynomial values for the respective fluxes. In this way, we naturally ended up with an ensemble LC that had no outliers, but showed some trace of the ‘trimming’ made. Figure 3 displays where the polynomial replacement of the original data points were made (green dots). When both the target and all of the comparison stars had to be corrected, we see a continuous sequence of points. In all other cases, the corrected points scatter around the ridge represented by the binned LC (yellow points).

The linear trend and the jump in the ensemble LC (issues 2 and 3) were treated within an iterative process by filtering out these systematic errors, fitting the cleaned LC to an eclipse model, and then subtracting this eclipse model from the starting dataset to get the next approximation for the systematic errors.

The systematic errors were represented by a linear function for the trend and a jump function to handle the discontinuity mentioned above: (1)

where F is the observed flux, t is the time measured from the first data point, and H is the Heaviside function with unit step at t_jump = 0.068 d. The jump position was fixed throughout the fit. Because the star blocks all radiation from the planet, the trapezoidal approximation for the occultation light curve suits perfectly: (2)

where t₁, t₂, t₃ and t₄ are the moment of ingress (first contact), the start and end of the total eclipse, and the moment of egress, respectively. The length of the ingress and egress phases are assumed to be equal: t₂ − t₁ = t₄ − t₃ = Δt. Except for the eclipse depth δ, all these parameters are scanned for the best fit within the framework of robust least squares. For any given set of {t_i} the transit depth was fitted in one step because of the linear nature of the parameter. The systematic error parameters {c_i} were fitted in the same manner. The final light curve for set-2 is shown in Sect. 3.2

Fig. 2 Raw flux variations for set-2 (see Table 1) and the resulting ensemble light curve (target flux over reference flux, ~ F1 /(F2 + F3), without outlier correction, normalized to its average). Stars 2 and 3 (see Fig. 1) serve as comparison stars. The eleventh-order polynomials, robustly fitted to the fluxes to handle outliers, are shown by green lines. The shaded area in the bottom panel indicates the period of intermittent clouds. The binned light curve (with overlapping bins – see text) is shown by yellow dots. The time axis is shifted to the moment of the first data point (BJD1).

Fig. 3 Ensemble LC obtained from set-2 after polynomial outlier correction. Black dots are original data points. Green dots are from the polynomial corrections. Yellow dots are from binning the data points independently of their origins (corrected or not corrected).

Fig. 4 Three light curves of Table 1 in the pre-merging phase. All light curves are filtered out from time-dependent linear trends and are normalized by the total (star+planet) flux Ft. Black dots are the binned values and yellow lines are the trapezoidal models fitted to the original (unbinned) data shown by deep gray dots. For better visibility, we increased the point size for set-1.

3.2 The Three Light Curves

In trying to treat all three datasets in the same way, that is, by starting from the simple ensemble light curve and employing the “minimum massage” post-processing step, we found that the case of set-1 (Chen et al. 2014) is different. As the ensemble light curve suffers excessively from systematic errors (see Fig. 2 of that paper), we decided to use their processed light curve that was obtained by applying carefully chosen external parameters (such as stellar position and image size) to separate systematic errors. On the other hand, for set-3 of Zhou et al. (2015) we used their simple ensemble light curve, even though there is also a substantial nonlinear trend in the data. In spite of this, we decided not to use an airmass or some polynomial correction (as given in the original paper), because this may introduce unpredictable changes in the eclipse and significantly depress the depth of the occultation.

All three light curves (sets-1 and -3 as above, set-2 as given in Fig. 3) serve as the input time series to fit them individually by the eclipse model and a linear trend (extended by a jump function for set-2). The result of this procedure is shown in Fig. 4.

Fig. 5 Trapezoidal models fitted to the Ks observations and folded by the orbital period. The shaded rectangle shows the expected event of occultation, assuming circular orbit and using the updated orbital period and transit center: P = 1.6284300 d, Tcen = 2458355.50805 [BJD] (see Sect. 4).

3.3 Merging the Three Light Curves

Before constructing the merged light curve, we need to check whether or not such a merging is possible, that is, if there is a unique orbital period that matches all three light curves within the observational errors. An evaluation of the updated orbital period using the primary transit observations from the TESS satellite and combination of the ephemerides with earlier followup data is presented in Sect. 4. Here we merely use the orbital period and the moment of the transit derived from that analysis.

By fitting the individual folded light curves, we can examine whether or not the data suggest strong discrepancies which would require consideration during the merging process. Figure 5 shows the individual phase-folded fits, indicating that the three datasets are in reasonable agreement, even if we consider set-1, the most discrepant of all. Set-1 contains the least number of data points (699 vs. 2084 and 1732 for sets 2 and 3, respectively), and also has the largest residual (data minus fit) scatter (in relative flux units: σ = 0.0036, vs. 0.0030 and 0.0033). In spite of these differences, all three datasets yield remarkably close egress phases. The reason for this is not entirely clear at this moment. In some cases, it might simply be the sign of more stable sky conditions in the second part of the run (which was indeed the case for set-2).

In the final step of the merging process, we packed all data points in a single phase-folded dataset by discarding the relatively small differences in data quality (i.e., weighting all data points from all sets equally)⁶. The phase-folded light curve, containing all the 4515 data points, was robustly fitted by the trapezoidal model. The resulting binned light curve and the best-fit trapezoidal are shown in Fig. 6. The fitted parameters are listed in Table 2. The errors were computed from simple Monte Carlo simulations, whereby the binned light curve (mapped back to all the 4515 phase points) was perturbed by a bin-dependent Gaussian noise. We opted to use the binned time series rather than the trapezoidal fit, because of the remaining systematic errors, especially before the ingress and after the egress. We generated 500 mock time series, fitted trapezoidals to each realization, and, after completion, we computed the standard deviations of the parameters. We refer to these standard deviations as the 1σ errors of the respective parameters. In Appendix A, we provide further details of the error calculation and the improvement of the parameters using the merged data as compared with the fits to the individual datasets.

Fig. 6 Upper panel: phase-folded occultation curve of all the Ks observations with bin averages (black dots). Lower panel: binned light curve of the above dataset with the errors of the bin averages and the trapezoidal model (black line) fitted to the original (unbinned) data shown in the upper panel.

4 Analysis of the TESS Data

We use the light curves acquired by the full sky survey satellite TESS for: (a) updating the ephemeris of the transit (as the occultation and the transit data were acquired in different epochs, we need a precise ephemeris to predict the transit phase immediately prior to the occultation if we want to estimate the eccentricity); and (b) measuring the emission in the optical, for which we need space-based data because of the high precision needed to detect the signature of the planet at this wavelength in the phase of occultation; the thermal emission in the optical is small because of the low temperature of the planet and the albedos of the gas giants in general are also small (e.g., Wong et al. 2020, 2021).

WASP-5 was observed by TESS in sector 02 between August 22 and September 20, 2018. The object was then revisited while scanning sector 29 between August 26 and September 22, 2020. The two segments comprise over 30 000 data points altogether in the short cadence (2 min) sampling rate. We note that Wong et al. (2020) previously performed an analysis of the sector 02 data and ended up with similar conclusions to ours, as we discuss in the subsections below.

Figure 7 shows the light curves from the above two sectors after employing the Presearch Data Conditioning (PDC) method of Smith et al. (2012) and Stumpe et al. (2012) implemented in the TESS pipeline⁷. The Simple Aperture Photometric (SAP) time series served as the input for the PDC filter. Both types of data were downloaded from the STScI MAST site⁸. We filtered the data further using a 36th-order robust polynomial fit to minimize the effect of the remaining systematic errors and possible stellar variability. The effect of this filtering is discussed in Appendix B.

Because sector 02 SAP data suffer from a large number of outlying data points, we are required to performed an iterative 3σ clipping for all datasets in order to make the analysis uniform. The clipping was made relative to the transit model and the clipped values were set equal to the corresponding model values.

Fig. 7 TESS light curves of WASP-5 from the sectors shown in the bottom right corners. The light curves have been processed by the TESS pipeline using PDC systematic-error corrections on the SAP fluxes. The 36th-order polynomial fit is shown with a green line, and is used to filter out the remaining systematic errors and possible stellar variability.

4.1 Updating Transit Ephemeris

To derive transit light curves free from other variations, we employed the same type of robust iterative method as briefly described in Sect. 3.1. The input data were the PDC/SAP time series as mentioned above. The model time series constitutes two multiplicative parts: the transit and a 36th-order polynomial. For the transit, we adopted the simple model of Kovacs (2020), representing the ingress and egress phases as linear flux depressions with the same steepness and duration. The limb darkening was modeled by a scalable U-shaped function. We found this model quite satisfactory at the level of the accuracy of the data analyzed.

The transit parameters for the various time series are shown in Table 3. After combining these with the transit parameters obtained from the five followup observations of Baluev et al. (2019), Moyano et al. (2017), Hoyer et al. (2012), Fukui et al. (2011), and Anderson et al. (2008), we found that the orbital period of Fukui et al. (2011) should be decreased by 0.123 s to properly match the published epochs⁹. By choosing sector 02 timing as a reference, the final ephemeris is given in Table 4. The error of the epoch was computed from 50 simple Monte Carlo simulations using the PDC data and is equal to the standard deviation of the epochs obtained from the 50 realizations. The error on the period was calculated from , where σ₁ = 0.00019 d as given in Fukui et al. (2011), σ₂ is the epoch error as given in Table 4, and the integer in the denominator is the elapsed epoch number between the two epochs. It is worth noting that the currently published ephemerides by Ivshina & Winn (2022) are in complete agreement with ours. There are 7 ms and 17 s differences between the periods and transit centers, respectively, corresponding an agreement within 1–2σ.

4.2 Search for Occultation and Phase Variation

To test the dependence of a possible detection of these delicate features on the data-processing methods, we used four data types: an SAP light curve (see Sect. 4.1) and a PDC light curve, both with and without robust polynomial correction. After prewhitening by the transit, we performed a bin signal search in the light curve folded by the orbital period. To account for the other possible variations, we employed a fully binned analysis, where the out-of-eclipse region was also divided into bins of the same size as the eclipse duration. After the bin with the largest flux depression was identified, we used a simple statistic to characterize its significance. Similarly, the phase variation was studied simply by a single-component Fourier fit to the original (i.e., not binned) phase-folded light curve. Significance tests were performed using injected signals into pure Gaussian time series. Further details on the secondary eclipse and phase variation searches together with the supplementary statistical tests are given in Appendix B. Here we summarize the constraints derived in that Appendix.

First, for the illustration of the data quality at the expected level of reflected light variation, we show the transit- and polynomial-filtered PDC-processed light curve in Fig. 8. The blue dots resulted from overlapping binning (see Sect. 3.1) with 200 bins (400 points altogether). The bin model has a bin width equal to the transit length, yielding 17 bins. Although this figure shows little indication of the presence of the type of signals we are searching for, as shown in Appendix B, the parameters fitted to the data of various processing levels remain remarkably stable. This leads to the following average values of the secondary eclipse depth and phase variation amplitude: δ = 0.157 ± 0.056 ppt, 2A = 0.113 ± 0.041 ppt.

With additional statistical tests, we found that at the observed amplitude of the cosine component, there is only a 0.3% probability that the underlying phase variation has a total (peak-to-peak) amplitude of greater than 0.20 ppt. Also, for a boxy eclipse of the same depth, the probability that the bin model yields a phase solution outside the expected secondary eclipse phase is less than 10%. With the observed correct location of the main dip for all four datasets, this suggests that we may have found a signature of the underlying signal.

Although the phase variation seems to yield a more stringent limit on the eclipse depth, the discordant phase¹⁰ of the cosine fit is discouraging, and we are reluctant to rely on the result suggested by this fit. Therefore, we use the eclipse depth quoted above as our best guess at present for the real secondary eclipse depth in the TESS waveband.

Fig. 8 Phase-folded and binned light curve of WASP-5 for the full TESS dataset after dividing the PDC flux values by the transit and polynomial components. The upper and lower panels, respectively, show the bin and cosine models (black lines). The transit phase is shifted to −0.5 for better visibility of the neighborhood of the expected phase of the secondary eclipse.

Fig. 9 Injected transit test of the full TESS dataset. We used the SAP light curve to inject the signal and then robust polynomial filtering was employed to lower the red noise. We show the frequency spectrum of the time series derived in this way, after subtracting the transit signal of planet b. Red arrow indicates the peak due to the injected signal with a transit depth of 0.3 ppt. The inset shows the close neighborhood of the test signal.

4.3 Search for Additional Transit Components

Although hot Jupiters systematically avoid close planetary companions (Poon et al. 2021), whether or not WASP-5 is one of those rare systems is still a matter of interest. Unfortunately, the short time-spans of the TESS observations complicate the search for the more common longer period companions, leading to lower observed multiple system rates from the TESS data (Otegi et al. 2022).

After prewhitening by the transit, we performed BLS searches (Kovacs et al. 2002) in the frequency interval [0.01, 10] c/d. The time series contains two dense tracks separated by ~730 days, comprising 36778 SAP data points altogether. We tested all four data-type combinations (SAP, PDC, both with and without polynomial filtering). All data types show an increasing power excess from 1 c/d down to 0.01 c/d with no prominent peak in this frequency interval. The spectra are flat in [1, 10] c/d, without any dominant peak superposed on the white noise background.

To test the detection limit in the potentially interesting frequency interval of low-order resonance, we injected a transit signal in the original SAP time series. We then performed a polynomial filtering as mentioned earlier in this section. The injected signal had a period of half of the orbital period of planet b and a transit depth of 0.3 ppt (corresponding to 1.8 Earth radii). We used a boxy transit with the same duration as that of planet b. The result is shown in Fig. 9.

From the structure of the spectrum, it is clear that 0.3 ppt transit depth is close to the lower limit of a transit signal we can hope to detect in the available dataset. This limit changes as a function of dataset and frequency, and is clearly higher for signals with periods longer than one day.

5 The Eccentricity

With the occultation ephemeris derived in Sect. 3 and with the transit ephemeris updated using the TESS data in Sect. 4, we can easily compute the two components of the eccentricity. For an easier reference, the components are as follows (Winn 2014): (3) (4)

where φ_obs and φ_cal are the observed (corrected for light-time effect) and calculated phases of the occultation centers (the latter is with the assumption of circular orbit), respectively. The argument of periastron is denoted by ω. As described in Sect. 3.3, the errors of the occultation signal in the Ks band were computed from a simple Monte Carlo simulation based on the binned version of the merged data from the three data sources. The observational noise was considered to be multiplicative and non-stationary, according to the standard deviations around the bin means. For the 500 realizations, we calculated the eccentricity components from the fitted trapezoidal occultation parameters. As is obvious from Eq. (3), the two components are not independent. Furthermore, e cos ω is expected to be less noisy, as several authors noted previously (e.g., Winn 2014). Indeed, Fig. 10 clearly shows both the correlation and the considerably tighter behavior of e cos ω. We also observe that the e sin ω component is shifted to more positive values. This is because noise makes the ingress/egress parts shallower, leading to the preference of longer eclipse durations in the best-fit search.

From these simulations, we obtained the errors also for the eccentricity components and, finally, for the eccentricity. For a simple reference, we summarize these parameters in Table 5. The correlation between the two eccentricity components is also exhibited by the lower error on the eccentricity than that on the e sin ω component.

It is instructive to compare our eccentricity values with those derived from the Spitzer data by Baskin et al. (2013). First we checked if there was any difference between using our transit ephemerides and those employed by Baskin et al. (2013). We found that for the two epochs Baskin et al. (2013) published in their Table 1, our ephemerides predicted an average offset 1 min greater than the one calculated from the ephemerides of Fukui et al. (2011) (3.7 min vs. 4.7 min). Their offset time implies e cos ω = 0.0025 ± 0.0012. Because the agreement is at the ~1σ level between their value and ours, we can average them out and arrive at a value of 0.0020 ± 0.0016, implying no difference from zero eccentricity¹¹.

Fig. 10 Eccentricity components obtained from the Monte Carlo simulations as described in the text. Gray and black dots, respectively, denote the e sin ω and e cos ω components. The inset shows the correlation between the two components, leading to a smaller error on the eccentricity as compared to that of the e sin ω component alone.

6 The Emission Spectrum

Here we examine how the more accurate occultation depth in the NIR (Sect. 3.3) and our preliminary estimate on the same quantity in the visible from the TESS data (Sect. 4.2) can constrain the atmospheric properties of WASP-5 b. The secondary eclipse analysis was presented in Sect. 3, where we derived an occultation depth of δ(occ, Ks) = (2.70 ± 0.14) ppt in the NIR.

In the visible, corresponding to the wide-band filter of TESS¹², we use the average of the eclipse depths obtained from four types of datasets: δ(occ, υis) = (0.157 ± 0.056) ppt. Although we could also use the value obtained from the estimation of the phase variation, we opted not to use this value for reasons discussed in Sect. 4.2.

Baskin et al. (2013) measured the emission at 3.6 µm and 4.5 µm from WASP-5 b using data from the Spitzer infrared satellite. Although we do not make any model fitting in this paper because we use the same models as given by Chen et al. (2014), we find it instructive to display all currently available data on the same plot.

The atmospheric models presented by Chen et al. (2014) are based on the plane-parallel equilibrium models of Madhusudhan & Seager (2009, 2010), employing free pressure–temperature profile and chemical composition. Figure 11 shows these theoretical spectra and the black body lines for fully efficient and zero circulations (i.e., lack of heat exchange between the day and night sides of the planet; see Cowan & Agol 2011). The atmospheric models have a monotonic pressure–temperature profile (i.e., no temperature inversion). The depth of the atmosphere was chosen to fit the brightness temperatures corresponding to the J, H, and K data of Chen et al. (2014).

There are two essential conclusions we can draw from the positions of the new data points in respect to these models. First, the lower error bar on Ks increased the significance of the higher Ks flux and suggests strong emission in this waveband. This emission could be due to additional emitters at a deeper level of the atmosphere (corresponding to temperatures higher than 2700 K).

Second, even though our estimate for the secondary eclipse depth is only tentative in the TESS band, it still yields a useful piece of information. This is because of the relatively small error on the data with respect to the spectral features in the optical. In particular, the current value of the emission in the visible corroborates what the Spitzer data may also indicate, that is, no strong preference for any of the models used.

From the optical occultation depth we can also estimate the geometric albedo. For an easier reference, here we repeat the necessary formulae presented by Cowan & Agol (2011) and, for example, by Daylan et al. (2021). The observed occultation depth constitutes two parts: the thermal radiation by the planet and the reflected light of the host star (5)

Assuming a circular orbit, the reflected light is directly related to the geometric albedo A_g by (6)

where R_p is the planet radius and a is the semi-major axis. The dayside thermal emission can be parameterized as (7)

where R_s is the stellar radius, and F_p and F_s are the wavelength(λ)-dependent fluxes of the planet and the star, respectively. In the black body approximation, parameter α is used to relate the substellar temperature to the dayside temperature T_day: (8)

The single parameter α comprises the Bond albedo A_b and the atmospheric circulation parameter ϵ. We refer to Burrows et al. (2008) and Cowan & Agol (2011) for the origin of the coefficients in the expression of α. Separation of A_b and ϵ is not possible using the occultation alone. However, by measuring the phase curve, one may attempt to derive the night side emissivity that depends solely on ϵ (Singh et al. 2022). Because of the lack of very high-quality data required by this method, we opted to parameterize the derived geometric albedo depending on the extreme limits of ϵ and omitting the negligible temperature decrease due to the expectedly small Bond albedo (e.g., Mallonn et al. 2019). Furthermore, as usual, we assume pure black body radiations both for the star and the planet in the waveband of interest.

The result is shown in Fig. 12. Although with fully efficient circulation (α = 1/4) the geometric albedo can be as high as A_g = 0.43 ± 0.15, based on the Ks occultation data (e.g., Kovács & Kovács 2019) and several other, more direct studies (i.e., those based on full phase curve analyses, such as Keating et al. 2019) it is highly unlikely that WASP-5 b stands out from the other hot Jupiters, that mostly have low circulation efficiency. Therefore, it is quite reasonable to assume that the true value of A_g is closer to the no circulation limit of A_g = 0.29 ± 0.15.

Fig. 11 Model thermal emission spectra from Chen et al. (2014) (gray for oxygen- and black for carbon-enhanced chemical compositions). Black body lines for fully efficient (α = 0.25) and inefficient (α = 2/3) circulations are shown by blue lines. The corresponding temperatures are given in the lower right corner. The optical and infrared occultation depths derived in this paper are shown with red dots, and the Spitzer data of Baskin et al. (2013) with yellow dots. The inset zooms in the optical waveband.

Fig. 12 Share of the thermal and reflected light in the total flux change during occultation with varying geometric albedo. The dots are for two extremes of planetary atmosphere dynamics with negligible Bond albedo: α = 2/3 (complete lack of circulation), and α = 1/4 (fully efficient circulation). The 1σ error of the observed value is indicated by the gray-shaded stripe. The vertical stripes show the resulting geometric albedos.

7 Conclusions

In this work, we deal with the secondary eclipse (occultation) light curve of the hot Jupiter WASP-5 b. Our goal is twofold: (i) we aim to derive an accurate occultation light curve in the 2MASS Ks band, and (ii) to use the latest TESS data to obtain the first estimate of the occultation depth in the optical. For goal (i) we used previously published Ks photometry by Chen et al. (2014) and Zhou et al. (2015) and combined these with the observations published in this paper for the first time. These observations were acquired by the FourStar infrared imager of the Baade 6.5 m telescope in 2011. Following the principle of “minimum data massage” we ended up with a high-precision occultation light curve. The relative flux depression (planet vs star) is 2.70 ± 0.14 ppt, which places WASP-5 b among the top few extrasolar planets with a Ks occultation light curve of such high relative precision.

We attempted to find the signature of secondary eclipse and phase variation in the visible by using the currently available TESS data from two visits. Based on the statistical tests presented in Appendix B, we find convincing pieces of evidence that the underlying signals are unlikely to have total variation greater than 0.20 ppt. As a result, we accepted the eclipse depth derived from the formal fit, that is (0.157 ± 0.056) ppt. Using these values, our main conclusions on the atmospheric properties of WASP-5 b are as follows.

• Simple black body radiation fails to reach the observed Ks emission at the level of 10σ. A similar statement, with somewhat lower significance of 4–7%, is also true for the band-averaged values of the adopted atmospheric models. Detailed atmospheric modeling with strong emission features in the Ks band is required to fit the high observed emission.

• The value derived for the emission in the TESS waveband shows no preference for any of the adopted models with oxygen or carbon enhancements. The observed emission value is ~2σ apart from both models.

• From the TESS eclipse depth we find that, depending on the circulation model, the geometric albedo A_g is likely to be in the range of 0.29–0.43. This places WASP-5b among the most reflective extrasolar planets but with the caveat that the detection of the secondary eclipse in the optical is preliminary in nature.

It is also worth mentioning that the ephemeris of the Ks occultation light curve further confirms the low (likely zero) eccentricity of the orbit, namely e = 0.005 ± 0.015. Furthermore, the TESS data do not suggest the presence of any additional transiting planet larger than ~2 Earth radii with a period of between 30 and 0.1 days.

The Ks waveband is within relatively easy reach for most of the ground-based telescopes with NIR capabilities. However, the timescale and the signal level mean that extrasolar occultation measurement are still challenging. These challenges are compounded by the local conditions, particularly in the form of the red noise component of the observed photometric time series. Perhaps the best way of handling red noise (if there are no other ways to filter it out), is to take multiple samplings. Careful combination of these samples will reduce both the white and red noise components. Due to the sparse sampling from the side of the available data points in the different wavebands, measurement accuracy is crucial for spectral retrieval. Although this spectral band is (or will be) available in various space missions (JWST now and ARIEL by the end of the decade – see Tinetti et al. 2018), the expected high demand (in particular for JWST) makes ground-based observations still very important in supplying high-quality data for more reliable extrasolar planet atmosphere modeling.

Acknowledgements

Constructive comments by the referee are appreciated, in particular those that have led to a deeper investigation of the significance of the secondary eclipse in the TESS band. We thank Nikku Madhusudhan for the valuable comments regarding the atmospheric modeling of WASP-5b. 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. I.D. was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subproject A03). G.C. acknowledges the support by National Natural Science Foundation of China (Grant No. 42075122, 12122308). Support from the National Research, Development and Innovation Office (grants K 129249 and NN 129075) is acknowledged.

Appendix A Trapezoidal fits to the nightly Ks data

In this Appendix we give further details of the method of the error estimation of the trapezoidal fit. This section also describes our closer assessment of the improvement resulting from merging the three datasets available in Ks color.

The individual nightly data presented by Chen et al. (2014), this paper, and Zhou et al. (2015) are fitted by the trapezoidal model of Sect. 3.1. The procedure followed for the individual datasets is the same as described in Sect. 3.3. The starting time series are the light curves that have already been cleaned of systematic errors and major outliers (see Fig. 4). To include the effect of changing noise level (both random and left-in systematics) throughout the night, the light curves were binned and the average bin values were used as the starting light curve for the simple Monte Carlo simulations. These simulations were performed with 500 noise realizations superposed on this “noiseless” light curve. The standard deviation of the superposed Gaussian noise changed according to the standard deviation of the observed data in each bin.

We used 80 overlapping bins, i.e., the total time-span of a given dataset was divided into 40 equal time segments and the averages and standard deviations of the data points belonging to each of these segments were calculated. The overlapping sequence started at half of the first bin and continued until the half of the last bin. The bin width of the second sequence was somewhat shorter, because we generated the same number of bins (i.e., 40) as for the first sequence. Finally, all the 80 bin averages and standard deviations were assigned to the 80 bins of equal length filling in the full time-span. This method of bin generation was used for all three datasets, independently of their lengths. While generating the mock light curves for noise estimation, we mapped back the bin statistics to the original timebase by choosing the bin values at the given time (or phase) value of the unbinned time series.

The simulated light curves served as inputs for our robust trapezoidal fitting routine. The resulting parameters were used to compute their standard deviations (1σ errors). The trapezoidal parameters themselves were obtained from a direct fit to the original (unbinned) data. In searching for the best fit, we used simple parameter scanning in the following ranges (in the units of the orbital phase): 0.45151 < T_ing < 0.49151, 0.04747 < T₁₄ < 0.07120, 0.00409 < T₁₂ < 0.00613, corresponding to ±0.02 phase change for T_ing with respect to a finally accepted best-fit value. The scanning ranges in T₁₄ and T₁₂ correspond to ±20 % relative variations (again, with respect to the same final best-fit values). Except for T₁₂, we had very rare hits at the limiting values. The more frequent hits for T₁₂ can be tolerated because of the limits imposed by the known system parameters.

The result is shown in Table A.1. To make the comparison easy, we copied the result of the analysis of the merged data from Table 2 in the fourth row. In the last but one column we show the standard deviation of the fit to the original data. The errors of the standard deviation of σ_fit come from the Monte Carlo simulations and indicate the data point number dependence of the statistical stability of the realizations.

We see that the merged data yield a considerable improvement with respect to all datasets. With an overall error decrease of ~ 20%, a relative modest improvement can be observed for set-2. For the other sets, we find improvement in the range of ~ 50%–300%. It is also observable that while our transit depth error is in agreement with the one derived by Zhou et al. (2015), there is nearly a factor of two difference between our error of 0.037 ppt and that of Chen et al. (2014). It is quite likely that the latter estimate (0.062 ppt) is closer to the real error, because their estimate follows the full sequence of light curve evaluation (which is certainly very sensitive to the method of systematics treatment), whereas our estimate relies on their “final product’’. However, from the point of view of the merged data analysis, the weight of this dataset is relatively small because of the small number of data points.

Appendix B Testing occultation and phase variation limits from TESS

We performed eclipse and phase variation analyses on the full set of TESS light curves – including sectors 02 and 29. With the option of polynomial filtering (see Sect. 4), we have four types of datasets to be analyzed (SAP and PDC, with or without polynomial filtering). The time series models were incomplete, i.e., they contained either an occultation or a phase variation signal. This is a reasonable approximation in the high-noise regime, when focusing on the detection limits of signals with very different time dependencies.

For the eclipse dip search, the folded light curve was approximated by a bin model, with equal width for all bins equal to that of the transit. The distribution of the bins was fixed, as given by the predicted phase of the secondary eclipse. The bin average yielding the largest flux depression and the center phase of that bin yields, respectively, the estimated eclipse depth δ and eclipse phase ψ. To characterize the quality of this representation of the eclipse, we used the dip significance parameter (DSP) as defined in Kovacs (2020) (B.1)

where e_δ is the error of the eclipse depth δ, σ_ab is the standard deviation of the bin averages for the bins outside the eclipse, σ_db is the standard deviation of the successive differences between the same bin averages. The last two terms are devoted to accounting for regular, smooth variations (σ_ab) and irregularities originating from pure noise at the phase scale of the bins (σ_db). The DSP value associated with the deepest bin average (i.e., the largest δ) was used to characterize the quality of the eclipse solution for a given dataset. It is noted that for a complete signal, DSP is limited by the second term in the denominator, because of the near equality of the eclipse depth and the total amplitude of the phase variation. Figure B.1 shows the relation between the light curve quality and the above parametrization.

The phase variation search was performed by a simple two-parameter (A and ϕ₀) least squares fit of a cosine function A cos(2π(ϕ − ϕ₀)). With a good approximation, the amplitude A of the cosine is half of the occultation depth (e.g., Daylan et al. 2021). The quality of the cosine model is characterized by the signal-to-noise ratio S/N = A/σ_fit, which, as expected for low-S/N signals such as those studied here, is tightly correlated with A.

The results of the above separate fits for the four types of datasets are displayed in Table B.1. Although there are differences among the datasets, the properties of the signals are largely stable. The box fit yields remarkable coincidence with the expected occultation center. On the other hand, the maximum of the cosine is systematically shifted, and precedes the occultation center in phase by ~ 0.36. The amplitude of the phase variation and the eclipse depth are consistent within the error limit: 〈δ〉 = 0.157 ± 0.056 ppt vs 〈2A〉 = 0.113 ± 0.041 ppt.

Fig. B.1 Parametrization of the phase-folded, binned test signals by DSP (see Eq. B.1). Several realizations of the signal described in the text and in the caption of Fig. B.2 were used with injected box depth of δ = 0.15 ppt. The width of the bins is equal to the eclipse duration of WASP-5, and are distributed according to the eclipse phase (centered arbitrarily at phase zero). In both cases, one realization is shown by a black line to follow bin-by-bin fluctuations more easily. To avoid unnecessary jamming, we show only realizations with the largest negative bin average at the eclipse center. The difference in the overall eclipse depth in the two cases is due to the different DSP cutoffs used.

It is important to address the statistical significance of these near-noise-level detections. To do so, we performed simple Monte Carlo simulations by generating the following type of signals: (B.2)

where G(i) is an uncorrelated Gaussian noise with standard deviations shown in Table B.1. Function S(i) is either a simple box or a cosine function, representing the secondary eclipse or the phase variation. The box model is zero out of the eclipse and −δ within the eclipse, with the box centered at the phase of the assumed eclipse, φ₀ (0.5 phase after the transit). Probability density functions (PDFs) and cumulative distribution functions (CDFs) were computed for the phase differences (i.e., Δφ = ψ − φ₀), DSP, and A values.

The primary goal of this investigation is to pose an upper limit on the underlying signal using the near-oise-level values derived above. Therefore, we computed two basic quantities. In the presence of injected signals of varying amplitudes and box depths we computed: (a) the occurrence rates of box depths and cosine amplitudes less than the observed values; (b) the occurrence rates of the dip center and cosine maximum phases within the proximity of the predicted secondary eclipse.

As an example, for three different injected box signals, Fig. B.2 shows the PDF of the phase differences and the CDFs of the DSP values. The phases seem to converge quite quickly (i.e., even for boxes as shallow as 0.1 ppt, more than half of the cases hit the near proximity of the expected phase). It is interesting to note that there is surplus of occurrences at the edges in the phase distribution. It is especially visible in the pure noise (δ = 0.0) simulations. The reason of this excess is the smaller number of data points in the bins at the end-phases of the folded time series (the bin occupancy depends on the positioning of the occultation and the width of the eclipse). This leads to more fluctuating values at the edges, and therefore a higher chance of being selected as the “best-fit” for the box model.

The DSP values are less sensitive to a small underlying signal; they are still near the pure-noise values and only become more distinctive from these if the box depths become deeper than 0.15 − 0.20 ppt.

A similar test performed on the same type of dataset with injected cosine signals indicates the opposite effect (see Fig. B.3). The phase settles at a far lower pace but the amplitude of the signal becomes more quickly detectable. For example, if the underlying cosine signal had an amplitude of 0.1 ppt (corresponding to an eclipse depth of 0.2 ppt), then the probability that we can detect a component as low as observed is less than 0.003.

Fig. B.2 Injected box test of pure Gaussian noise with σ = 2.32 ppt and n = 33256 data points. Right column: CDF of DSP for various injected box depths δ. Vertical bars show the observed value (6th column of Table B.1 in row 4). Left column: PDF of the phase difference between the test box center and the calculated secondary eclipse center from the bin model for the same box depths as shown in the corresponding panel on the right. The bin width is the total transit duration, i.e., in phase units ϵ = 0.030.

We can assess the likelihood of the various underlying signals in the currently available TESS data using the data settings (i.e., data point distribution, N and σ) for the two extreme data types shown in Table B.1. These settings correspond to the PDC and SAP fluxes, respectively, with and without polynomial detrending. We use the detection power of the phase of the bin search and the amplitude sensitivity of the Fourier fit. The upper panel of Fig. B.4 shows the occurrence rate of the observed total amplitude as a function of the underlying (injected) signal amplitude. The lower panel exhibits the likelihood that the best-fitting bin phase is not in close proximity to the predicted phase.

These plots suggest that the underlying phase variation should have a total amplitude of less than 0.20 ppt with a probability of more than 99.7%, because otherwise there would be a probability of less than 0.3% of obtaining a cosine amplitude as small as that given in Table B.1. Although at a somewhat lower level of significance, this result is corroborated by the frequency of the correct phase hits in the box test. Depending somewhat on the data type, the probability of not hitting the correct phase for an underlying boxy eclipse signal of depth 0.20 ppt is between 4% and 9%. We note that the average of the eclipse depths derived from the observed data and the associated 1σ formal error yield an upper limit of 0.157 + 0.056 = 0.213 ppt, close to the high-significance limit obtained above. This, together with the stability of the observed eclipse phases, lends further support to the tentative detection of the secondary eclipse for WASP-5b from the TESS data.

Fig. B.3 As in Fig. B.2 but for the injected cosine test. The amplitudes a of the injected cosine functions were chosen to be compatible with the corresponding box test (i.e., 2a = δ).

Fig. B.4 Upper panel: Occurrence rate of the detection of cosine amplitudes δ/2 less than δ0/2 in the presence of injected cosine amplitudes dinj/2. The Gaussian components of the mock signals were generated using the standard deviations of the SAP and PDC time series (first and last rows in Table B.1). The δ0 values refer to the respective total amplitudes in the same table (i.e., 0.086,0.106 ppt for the SAP and PDC data, respectively). Lower panel: Occurrence rate of phase difference Δφ for the same noise models as above but injected by box signals of depths dinj and widths of 2ϵ = 0.060.

References

All Tables

All Figures

	Fig. 1 Field of the first chip of the FourStar infrared mosaic imager of the Baade telescope. North is to the left, east is to the bottom. Chips 2, 3, and 4 are located clockwise starting at the bottom of chip 1. The image size is 5.5′ × 5.5′. We used comparison stars Nos. 2 and 3 only.
In the text

Fig. 2 Raw flux variations for set-2 (see Table 1) and the resulting ensemble light curve (target flux over reference flux, ~ F1 /(F2 + F3), without outlier correction, normalized to its average). Stars 2 and 3 (see Fig. 1) serve as comparison stars. The eleventh-order polynomials, robustly fitted to the fluxes to handle outliers, are shown by green lines. The shaded area in the bottom panel indicates the period of intermittent clouds. The binned light curve (with overlapping bins – see text) is shown by yellow dots. The time axis is shifted to the moment of the first data point (BJD1).

In the text

	Fig. 3 Ensemble LC obtained from set-2 after polynomial outlier correction. Black dots are original data points. Green dots are from the polynomial corrections. Yellow dots are from binning the data points independently of their origins (corrected or not corrected).
In the text

	Fig. 4 Three light curves of Table 1 in the pre-merging phase. All light curves are filtered out from time-dependent linear trends and are normalized by the total (star+planet) flux Ft. Black dots are the binned values and yellow lines are the trapezoidal models fitted to the original (unbinned) data shown by deep gray dots. For better visibility, we increased the point size for set-1.
In the text

	Fig. 5 Trapezoidal models fitted to the Ks observations and folded by the orbital period. The shaded rectangle shows the expected event of occultation, assuming circular orbit and using the updated orbital period and transit center: P = 1.6284300 d, Tcen = 2458355.50805 [BJD] (see Sect. 4).
In the text

	Fig. 6 Upper panel: phase-folded occultation curve of all the Ks observations with bin averages (black dots). Lower panel: binned light curve of the above dataset with the errors of the bin averages and the trapezoidal model (black line) fitted to the original (unbinned) data shown in the upper panel.
In the text

	Fig. 7 TESS light curves of WASP-5 from the sectors shown in the bottom right corners. The light curves have been processed by the TESS pipeline using PDC systematic-error corrections on the SAP fluxes. The 36th-order polynomial fit is shown with a green line, and is used to filter out the remaining systematic errors and possible stellar variability.
In the text

	Fig. 8 Phase-folded and binned light curve of WASP-5 for the full TESS dataset after dividing the PDC flux values by the transit and polynomial components. The upper and lower panels, respectively, show the bin and cosine models (black lines). The transit phase is shifted to −0.5 for better visibility of the neighborhood of the expected phase of the secondary eclipse.
In the text

Fig. 9 Injected transit test of the full TESS dataset. We used the SAP light curve to inject the signal and then robust polynomial filtering was employed to lower the red noise. We show the frequency spectrum of the time series derived in this way, after subtracting the transit signal of planet b. Red arrow indicates the peak due to the injected signal with a transit depth of 0.3 ppt. The inset shows the close neighborhood of the test signal.

In the text

	Fig. 10 Eccentricity components obtained from the Monte Carlo simulations as described in the text. Gray and black dots, respectively, denote the e sin ω and e cos ω components. The inset shows the correlation between the two components, leading to a smaller error on the eccentricity as compared to that of the e sin ω component alone.
In the text

Fig. 11 Model thermal emission spectra from Chen et al. (2014) (gray for oxygen- and black for carbon-enhanced chemical compositions). Black body lines for fully efficient (α = 0.25) and inefficient (α = 2/3) circulations are shown by blue lines. The corresponding temperatures are given in the lower right corner. The optical and infrared occultation depths derived in this paper are shown with red dots, and the Spitzer data of Baskin et al. (2013) with yellow dots. The inset zooms in the optical waveband.

In the text

Fig. 12 Share of the thermal and reflected light in the total flux change during occultation with varying geometric albedo. The dots are for two extremes of planetary atmosphere dynamics with negligible Bond albedo: α = 2/3 (complete lack of circulation), and α = 1/4 (fully efficient circulation). The 1σ error of the observed value is indicated by the gray-shaded stripe. The vertical stripes show the resulting geometric albedos.

In the text

Fig. B.1 Parametrization of the phase-folded, binned test signals by DSP (see Eq. B.1). Several realizations of the signal described in the text and in the caption of Fig. B.2 were used with injected box depth of δ = 0.15 ppt. The width of the bins is equal to the eclipse duration of WASP-5, and are distributed according to the eclipse phase (centered arbitrarily at phase zero). In both cases, one realization is shown by a black line to follow bin-by-bin fluctuations more easily. To avoid unnecessary jamming, we show only realizations with the largest negative bin average at the eclipse center. The difference in the overall eclipse depth in the two cases is due to the different DSP cutoffs used.

In the text

Fig. B.2 Injected box test of pure Gaussian noise with σ = 2.32 ppt and n = 33256 data points. Right column: CDF of DSP for various injected box depths δ. Vertical bars show the observed value (6th column of Table B.1 in row 4). Left column: PDF of the phase difference between the test box center and the calculated secondary eclipse center from the bin model for the same box depths as shown in the corresponding panel on the right. The bin width is the total transit duration, i.e., in phase units ϵ = 0.030.

In the text

	Fig. B.3 As in Fig. B.2 but for the injected cosine test. The amplitudes a of the injected cosine functions were chosen to be compatible with the corresponding box test (i.e., 2a = δ).
In the text

Fig. B.4 Upper panel: Occurrence rate of the detection of cosine amplitudes δ/2 less than δ0/2 in the presence of injected cosine amplitudes dinj/2. The Gaussian components of the mock signals were generated using the standard deviations of the SAP and PDC time series (first and last rows in Table B.1). The δ0 values refer to the respective total amplitudes in the same table (i.e., 0.086,0.106 ppt for the SAP and PDC data, respectively). Lower panel: Occurrence rate of phase difference Δφ for the same noise models as above but injected by box signals of depths dinj and widths of 2ϵ = 0.060.

In the text