© C. Jiang et al. 2022

1 Introduction

Transit spectroscopy is currently the most accessible and most widely used technique for revealing the absorption and scattering features of exoplanet atmospheres. Based on transmission spectra, we can identify the atomic and molecular species, along with the properties of clouds or hazes, in planetary atmospheres (Seager & Sasselov 2000; Brown 2001; Madhusudhan 2019). When transit spectroscopy is carried out with ground-based telescopes, special attention must be paid to the systematics of telluric origins. For high-dispersion spectroscopy, the telluric absorption lines can be directly resolved and corrected (Redfield et al. 2008; Snellen et al. 2008, 2010). For low-dispersion spectroscopy, the long-slit or the multi-object observational mode can be used to simultaneously monitor the fluxes of the target and the reference star(s), and then differential spectrophotometry can be applied to reduce the telluric effects. These methods have been widely employed in ground-based transit spectroscopic surveys (e.g., Sing et al. 2012; Nikolov et al. 2016, 2018; Sedaghati et al. 2016, 2017; Chen et al. 2017, 2018; Espinoza et al. 2019; Kirk et al. 2019; Carter et al. 2020).

We used spectroscopic light curves from ground-based observations to validate the methods commonly used in low-dispersion transmission spectroscopy. Our aim is to investigate to which extent the light-curve systematics would impact the derived transmission spectra. We performed this benchmark test by observing a white dwarf transiting a main-sequence star. The geometric features of transit light curves for white dwarfs are similar to those for hot Jupiters, except for the additional systematic effects from eclipsing binaries, including ellipsoidal light variation, mutual illumination, and Doppler boosting (Carter et al. 2011). According to Di Stefano (2011), the hydrogen envelop of the white dwarf progenitor is stripped off in the late phase of mass transfer with its main-sequence companion, leaving a degenerate helium core that is the white dwarf. Due to the extremely high surface gravities of white dwarfs, their atmospheric scale heights are so small that the transmission signals are undetectable in the visible to near-infrared wavebands, limited by current instrumental photometric precision, thus producing featureless transmission spectra.

We selected two transiting white dwarfs as our targets of interest: KIC 10657664B (also known as KOI-964B), and KIC 9164561B. They both transit A-type dwarf stars, and their transit parameters are close to those of typical hot Jupiters. The binary system KIC 10657664 was selected for its well-determined transit parameters, where the white dwarf companion orbits the A-type host star (2.23 ± 0.12 M_⊙, 1.89 ± 0.14 R_⊙, 13.64 R-mag) at a distance of about seven times the host star radius, with a period of ~3.27 days (Carter et al. 2011; Wong et al. 2020). Several reference stars were aligned, together with the target, inside the long slit of the spectrograph, which enabled us to investigate whether the selection of reference stars would impact the calculated transit parameters and transmission spectra. The other target, KIC 9164561B, orbits the primary star (2.02 ± 0.06 M_⊙, 2.54 ± 0.03 R_⊙, 14.06 R-mag) at a closer distance of ~2.5 times the host star radius with a shorter period of ~1.27 days (Rappaport et al. 2015). This binary system was selected for its oscillating light curves that arise from the stellar pulsation of the A-type host, which can be treated as an analog to other transiting exoplanets whose host stars also exhibit pulsation features in the light curves (e.g., WASP-33b; von Essen et al. 2014, von Essen et al. 2019). Based on this target, we investigated whether an unbiased transmission spectrum can be derived when the light curves are contaminated by considerable stellar pulsations.

This paper is organized as follows. In the next section, we summarize the transit observations and the data reduction procedures. We introduce the methods of light-curve analyses in Sect. 3. The results of benchmark tests are shown in Sect. 4. We draw our conclusions in Sect. 5.

2 Observations and Data Reduction

The transit events of our two targets were observed using the long-slit mode of the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy at the Gran Telescopio Canarias (GTC OSIRIS) with the R1000R grism. Two red optimized 2048 × 4096 Marconi CCDs were used to record images with a gap of 9.4″ in between. The pixel scale of the imaging was 0.254″ after a 2 × 2 pixel binning. The gap between the two CCDs was parallel to the dispersion direction and did not affect the OSIRIS spectra. The unvignetted field of view for imaging was 7.8′ × 7.8′. The slit was 7′ long and 12″ wide. The total waveband of the R1000R grism is 5100–10000 Å with an instrumental resolution of 2.6 Å per pixel. The standard 200 kHz readout mode was adopted, and the corresponding readout noise was ~4.5e⁻. The arc lines of HgAr, Xe, and Ne were measured through a 1″ wide slit with the same R1000R grism for the wavelength calibration of the science spectra.

The transit of KIC 10657664B was observed on 16 July 2015. During the observation, the Moon was absent and the weather was mostly clear. Four neighboring reference stars (hereafter REF-A1, REF-A2, REF-A3, and REF-A4) with similar magnitudes were observed simultaneously with the target star. Table 1 lists the identifiers, stellar types, and R-band magnitudes of all observed stars. The target and the REF-A1 star were placed on CCD1, while the other three stars (REF-A2, REF-A3, and REF-A4) were placed on CCD2 (Fig. 1).

The transit of KIC 9164561B was observed on 10 July 2015. The weather was clear during the observation, while the Moon was illuminated by ~36% at a distance of 87°. One reference star KIC 9099798 (hereafter REF-B1) was simultaneously observed with the target star, and both of them were placed on CCD1 (see Fig. 1). Although REF-A3 and REF-B1 are rotationally variable stars, their variation periods (~ 14.6 days for the former and ~13.2 days for the latter; McQuillan et al. 2014) are much longer than the duration of the observations, thus hardly affecting the differential flux correction. More details of the two observations are listed in Table 2.

The spectral data were reduced using the standard IRAF routines together with our customized IDL scripts based on AstroLib¹. We first performed corrections on overscan, bias, and flat field for all spectral images. We then constructed a two-dimensional pixel-to-wavelength mapping based on the row-by-row line identification of the arc lines. The science spectral images were transformed to the wavelength space so that the sky emission lines are no longer curved in the spatial direction. We masked out all the stars and fit the remaining sky backgrounds with a linear model to predict the variation in sky emission within the stellar point spread function (PSF). Then we transformed the 2D sky models back to the pixel space before subtracting them from the science spectral images. In addition, we performed a simple sigma-clipping method on the time series of each pixel to flag the points that were hit by cosmic rays, which were then replaced by the median values of the neighboring pixels in the corresponding exposures. Finally, we extracted the stellar spectra of the target and the reference stars with the optimal extraction algorithm (Horne 1986) through the aperture diameters ranging from 4 to 44 pixels (1.02″–11.18″) with an increment of 1 pixel. Since the target and the reference stars were not perfectly aligned to the central line of the slit and the stars drifted within 1 pixel during the observations, we first applied wavelength correction for each object and each exposure by calculating the wavelength shifts of the Hα and Na I lines from the laboratory air wavelengths. We then aligned the telluric absorption lines of the reference star to those of the target star based on the cross-correlation function of spectral profiles around the oxygen A-band (759–770 nm).

Although the R1000R grism was used in combination with a spectral order sorter filter, which cuts out the light blueward from ~495nm², there is a slight second-order contamination in the stellar spectra at 960–980 nm, especially for the two science targets with strong blue radiation. In addition, the fringing in the OSIRIS CCD images are measured to be ~5% for wavelengths larger than 930 nm. Therefore we only used the wavelength range of 515–915 nm for the light-curve analyses (Fig. 2). Integrating the spectra in the broadband (excluding 755–775 nm where the telluric oxygen A-band absorption is predominant), we constructed the raw light curves of all the stars. The other telluric absorption bands, including the oxygen B-band (~687 nm) and the water bands (e.g., ~820 nm), were not excluded because they do not have significant impacts on the resulting broadband light curves. The spectroscopic light curves were constructed by integrating the fluxes in uniform 20 nm passbands (Fig. 2), which are compromised between the signal-to-noise ratio (S/N) and the spectral resolution. We reduced the telluric effects by dividing the light curves of the target star with that of a reference star. The transit light curves were then normalized based on the out-of-transit flux. We determined the best aperture size by minimizing the light-curve scatter, which varies with different pairs of the target and the reference stars. For KIC 10657664, the best aperture diameters are found to be 10 pixels (2.54″) when using the REF-A1, REF-A3, and REF-A4 stars, and 26 pixels (6.60″) when using the REF-A2 star. For KIC 9164561, the best aperture diameter is 23 pixels (5.84″).

Fig. 1 Through-slit images of the two observations. Upper row: Positions of KIC 10657664 and the reference stars (REF-A1 to A4) in CCD1 and CCD2. Lower row: Positions of KIC 9164561 and the reference star (REF-B1) in CCD1. The red mark is the GTC Nasmyth rotator center. Only a part of the data section is displayed for clarity.

Fig. 2 Time-averaged stellar spectra of all objects. The upper panel shows the spectra of KIC 10657664 and its four reference stars. The lower panel shows the spectra of KIC 9164561 and its single reference star. The green shaded intervals indicate the wavebands for spectroscopic light curves.

3 Transit Light-Curve Analysis

3.1 Transit Model

The modeling of transit light curves is similar to our previous work (Jiang et al. 2021). We used the Python code batman (Kreidberg 2015) to calculate the transit models proposed by Mandel & Agol (2002). We adopted the transit ephemeris of Wong et al. (2020) for KIC 10657664 and those of Rappaport et al. (2015) for KIC 9164561 to roughly predict the central transit time assuming circular orbits. We used the quadratic stellar limb-darkening law in the transit model. The limb-darkening coefficients u₁ and u₂ were constrained by Gaussian priors, where the prior means were calculated with the code developed by Espinoza & Jordan (2015) based on the ATLAS9 model atmospheres (Castelli & Kurucz 2003), and the prior uncertainties were assumed to be 0.1. For broadband light-curve modeling, the free parameters in a transit model include the companion-to-host radius ratio R₂/R₁, the quadratic limb-darkening coefficients u₁ and u₂, the semimajor axis relative to the host star radius a/R₁, the orbital inclination i, and the central transit time t_c. The spectroscopic light curves were separately fit in each narrow band. Since the orbital parameters (a/R₁, i, and t_c) are independent of wavelength, we fixed these parameters to the median estimates derived from white-light curve fitting. Thus the free parameters used for spectroscopic light-curve fitting are R₂/R₁, u₁, and u₂.

The nightside thermal emission of most exoplanets is negligible with current instrumental precision, except for some ultrahot Jupiters observed in infrared wavebands (e.g., Kipping & Tinetti 2010; Chakrabarty & Sengupta 2020). In contrast, the self-emission from white dwarfs can significantly dilute the transit signals in the optical waveband, which would cause an underestimate of the transit depth if not corrected. The diluted transit model is calculated to be $$m^{\ast }\left(t;ℱ\right)=\frac{m\left(t\right)+ℱ}{1+ℱ},$$(1)

where m(t) is the transit model without the dilution effect, and ℱ is the companion-to-host flux ratio in the corresponding waveband. The flux ratio ℱ can be directly measured through secondary eclipse observations, but we lack such spectroscopic measurements covering the same wavebands. The flux ratio varies with wavelength and is degenerate with the transit depth, for which we cannot take the estimates from the Kepler data, nor can we directly retrieve it from the GTC spectroscopic light curves. Therefore we estimated the spectroscopic flux ratios using the blackbody assumption, $$ℱ\left(R_2/R_1,T_1,T_2;{\lambda }_a,{\lambda }_b\right)={\left(\frac{R_2}{R_1}\right)}^2\cdot \frac{{\displaystyle {\int }_{{\lambda }_a}^{{\lambda }_{\text{b}}}B\left(\lambda ,T_2\right)\text{d}\lambda }}{{\displaystyle {\int }_{{\lambda }_a}^{{\lambda }_{\text{b}}}B\left(\lambda ,T_1\right)\text{d}\lambda }},$$(2)

where R₂/R₁ is the radius ratio, T is the effective temperature, the subscripts 1 and 2 refer to the primary and the secondary star, respectively, (λ_a, λ_b) is the wavelength interval, and B(λ, T) is the Planck function. We adopted the temperature estimates from Wong et al. (2020) and Rappaport et al. (2015), which are $T_1={9940}_{-230}^{+260}$ K and T₂ = 15 080 ± 400 K for the KIC 10657664 system, and T₁ = 7870 ± 150K and T₂ = 10410 ± 200K for the KIC 9164561 system. We find that the measurement uncertainties of T₁ and T₂ have negligible impacts on the light-curve fitting. Therefore we used fixed values instead of assuming Gaussian priors on T₁ and T₂ so as to reduce the number of free parameters in the transit model.

3.2 Baseline Variation Effects in Binary Systems

The effects of Doppler boosting (DB), mutual illumination (ILL), and ellipsoidal light variation (ELV) may introduce additional periodic features in the phase curve of a binary system, and thus might change the baseline of transit light curves. Both Wong et al. (2020) and Rappaport et al. (2015) used third-order harmonic series to account for these effects when analyzing the Kepler light curves. As introduced in Carter et al. (2011), the features of DB, ILL, and ELV mainly correspond to the sin(ϕ), cos(ϕ), and cos(2ϕ) terms, respectively, where ϕ(t) = 2π(t − t_c)/P is the orbital phase and P is the orbital period. Based on the results of Wong et al. (2020) and Rappaport et al. (2015), we can roughly evaluate the baseline variation amplitudes (BVA) induced by these systematic effects during a five-hour transit observation with the GTC OSIRIS instrument. The total astrophysical induced BVA is ~40 ppm for KIC 10657664, which is negligible compared to the photometric precision of ~200ppm. Therefore we did not consider these astrophysical effects when fitting the transit light curves of KIC 10657664. For the other system, KIC 9164561, although the BVA due to the DB effect is only ~10 ppm, those due to the ILL and ELV effects can reach considerable amplitudes of ~300 ppm and ~4000 ppm in the GTC transit observation. This is because the amplitudes of ILL obeys the inverse square law and that of ELV is inversely proportional to the fourth power of the relative orbital distance. Since the value of a/R₁ for KIC 9164561 (a/R₁ ≈ 2.5) is considerably lower than that for KIC 10657664 (a/R₁ ≈ 7.0), the effects of ILL and ELV are consequently much stronger for the latter. Thus the significant baseline variations should not be ignored for KIC 9164561.

We corrected for the effects of ILL and ELV in the light-curve analyses of KIC 9164561 by multiplying its transit model with a harmonic series, $$h\left(\cancel{o}\right)=1+A_{\text{ILL}}\cos \left(\cancel{o}\right)+A_{\text{ELV}}\cos \left(2\cancel{o}\right),$$(3)

where the harmonic amplitudes A_ILL and A_ELV are free parameters. The other harmonic terms with small amplitudes, such as A_DB sin(ϕ), A_ELV cos(ϕ), and A_ILL cos(2ϕ), are not included in Eq. (3) because they may be treated as noise and be accounted by the noise models illustrated in Sect. 3.3. We took $A_{\text{ILL}}~\mathcal{N}\left(2385,36\right)$ (in unit of ppm) from Rappaport et al. (2015) for all the broadband and spectroscopic light curves of KIC 9164561. We do not expect a significant wavelength dependence of A_ILL according to the ILL model presented in Carter et al. (2011). We note that the ELV model depends on the host star limb-darkening and gravity-darkening effects (Morris 1985; Carter et al. 2011). Thus the wavelength dependency of A_ELV should be considered when fitting spectroscopic light curves. According to Eqs. (2) and (3) in Morris (1985), the amplitude of ELV is predicted to be $$A_{\text{ELV}}\approx -\frac{45+3u}{20\left(3-u\right)}\left(1+\tau \right)q{\left(R_1/a\right)}^3{\sin }^{2}i,$$(4)

where u is the linear limb-darkening coefficient, τ is the gravity-darkening coefficient, q is the companion-to-host mass ratio, R₁/a is the host star radius relative to the orbital semimajor axis, and i is the orbital inclination. Taking q = 0.097 ± 0.004, R₁/a = 0.397 ± 0.003, and i = 71.59 ± 0.22 deg from Rappaport et al. (2015), we need to calculate u and τ in each waveband so as to obtain the prior constraints on A_ELV. The linear limb-darkening coefficients were calculated with the code of Espinoza & Jordán (2015) based on the stellar ATLAS model (Kurucz 1979), which we previously used to calculate the quadratic limb-darkening coefficients for the transit model. The gravity-darkening coefficients were calculated with Eq. (10) in Morris (1985). The predicted values of A_ELV in each narrowband are listed in Table A.1, and their estimated uncertainties are ~400 ppm according to error propagation of Eq. (4). We calculated the weighted means of u and τ with the instrumental response curve of the OSIRIS R1000R grism as the weights to obtain u = 0.43, τ = 0.72, and A_elv = −8466 ± 400 ppm in the broadband, which is slightly lower than the measured value of A_ELV = −8880 ± 35 ppm from Rappaport et al. (2015) because the Kepler waveband covers shorter visible wavelengths starting from 420 nm.

3.3 Noise Model

There are typically two ways to handle light-curve systematics. One way is to use nonparametric stochastic models such as Gaussian process (GP) regression introduced in Gibson et al. (2012) and Foreman-Mackey et al. (2017). A GP model estimates the correlation between the input variables and the systematics by constructing a covariance matrix with a kernel function, which provides an adaptive and robust characterization of systematic noise. The other way is to use parametric deterministic models, where the systematics are accounted for by a set of baseline functions (BFs). There are various forms of baselines that can be defined, for example, a polynomial series of the state vectors or a harmonic series of time. Exponential baselines have been used to analyze the light curves of the Hubble Space Telescope (e.g., Berta et al. 2012; Kreidberg et al. 2014). The “best” baseline model needs to be selected from a large number of presumed models and depends on the specific observational conditions and instruments. Gibson (2014) tested these two methods using simulated transit light curves. He suggested that multiple techniques should always be used to test whether the inference is dependent on subjective choices made. Here, we focus on the performance of these methods on real data from ground-based observations and investigate the model dependence of the derived transmission spectra.

3.3.1 Using Gaussian Processes

We use GP regression to account for the correlated noise in transit light curves, which was implemented with the celerite code developed by Foreman-Mackey et al. (2017). The kernels supported by celerite are mixtures of exponentials, which enable fast modeling of one-dimensional GPs. A detailed introduction to the GP framework is available in Rasmussen & Williams (2006). Gibson et al. (2012, Gibson et al. 2013a,b) introduced and discussed the specific applications of GP modeling in transit light-curve analyses.

In a GP model, the observed light curve f follows a multi-variate Gaussian distribution, $$f\sim \mathcal{N}\left(m^{\ast }\left(t;\theta \right),K\left(t;\phi \right)\right),$$(5)

where m* is the mean function with a parameter vector θ, K is the covariance matrix with a parameter vector φ, and t is the time vector. The mean function m* is exactly the diluted transit model defined in Eq. (1). The covariance matrix K characterizes the residual noise r = f − m* consisting of correlated noise and white noise. The elements in K are determined by a covariance function (also called kernel function) k(τ_ij; φ), where τ_ij ≡ |t_i − t_j| is the time difference between two data points.

In this work, we test two well-defined kernels. One is the 3/2-order Matérn kernel (hereafter the M32 kernel), which takes the form $$\begin{array}{l}k\left(\tau ij;\phi \right)={\alpha }^2\left(1+\frac{\sqrt{3}\tau ij}{\ell }\right)\exp \left(-\frac{\sqrt{3}\tau ij}{\ell }\right),\\ \phi =\left(\alpha ,\ell \right),\end{array}$$(6)

where α and ℓ are the amplitude and length scale of the correlated noise. This kernel is approximated with the Matérn32 term in celerite. The other kernel function we test is a stochastically driven damped simple harmonic oscillator (hereafter the SHO kernel). It is computationally efficient for modeling stellar variations (Foreman-Mackey et al. 2017). The SHO kernel exhibits different properties depending on its quality factor Q. In the small Q limit (0 < Q ≪ 1/2), it is over damping and presents no oscillatory features, while in the large Q limit (Q ≫ 1), it becomes $$\begin{array}{l}k\left(\tau ij;\phi \right)=S_{0}Q{\omega }_{0}exp\left(-\frac{{\omega }_0{\tau }_{ij}}{2Q}\right)cos\left({\omega }_0{\tau }_{ij}\right),\\ \phi =\left(S_0,Q,{\omega }_0\right),\end{array}$$(7)

where S₀ controls the amplitude of the oscillator, and ω₀ is the characteristic frequency. This kernel can be used to characterize the asteroseismic oscillations (e.g., Grunblatt et al. 2017; Pereira et al. 2019).

To account for the white noise in light curves, we added a diagonal jitters term to the kernel function such that $$\begin{array}{l}{k}^{\prime }\left(\tau ij;{\phi }^{\prime }\right)={\epsilon }^2{\delta }_{ij}+k\left(\tau ij;\phi \right),\\ {\epsilon }^2={\epsilon }_{ii}^2+{\epsilon }_{\text{ex}}^2,\\ {\phi }^{\prime }=\left({\epsilon }_{\text{ex,}}{\phi }_1,\dots ,{\phi }_n\right),\end{array}$$(8)

where δ_ij is the Kronecker delta, ε is the total white noise, ε_ii is the flux uncertainty propagated from photon-dominated noise, and ε_ex is a factor to account for white-noise underestimation.

3.3.2 Using Parametric Baseline Functions

When using parametric baseline functions b(t; φ), the transit light curves are fit by $$f\left(t\right)=m^{\ast }\left(t;\theta \right)b\left(t;\phi \right)+\epsilon ,$$(9)

where φ is the parameter vector of the baseline, and ε is the total white noise (Eqs. (8)). Different from the GP models, usually many similar BF models can be enumerated, and model selection must be conducted to determine the best model. In this context, we compared the Bayesian evidence among all the BF models when fitting white-light curves and selected the one with the highest evidence as the best model that we used in the subsequent spectroscopic light curve analyses.

For KIC 10657664, the systematics are mainly instrumental and telluric origins. We define the polynomials of state vectors up to the third order as the BF models. We used six state vectors, which are the time t, the target position drift x along the dispersion direction and y along the spatial direction, the full width at half maximum (FWHM) of the target PSF in the spatial direction w, the instrumental rotation angle a, and the airmass z, to build a total of 4096 (= 4⁶) BF models. For instance, one baseline function can be written as $b\left(t\right)=c_0+{\displaystyle {\sum }_{n=1}^3{c}_{t,n}{t}^n+{\displaystyle {\sum }_{n=1}^2{c}_{w,n}{w}^n\left(t\right)}}$, and φ = (c₀, c_t,1, c_t,2, c_t,3, c_w,1, c_w,2). The state vectors were normalized by their maximum norms before being added into polynomials. Figure A.1 shows the time series of the original state vectors.

For KIC 9164561, the systematics are mainly due to stellar pulsations. We find that the polynomial baselines composed of state vectors are unsuitable for its light-curve fitting. Therefore we followed von Essen et al. (2019) and used the sum of sinusoidal functions to account for the oscillatory systematics. The sinusoidal BF for the light curves of KIC 9164561 has the form $$b\left(t\right)=c_0+{\displaystyle \sum _{i=1}^m{A}_i\sin \left(2\pi t/T_i-{\cancel{o}}_i\right),}$$(10)

where A_i is the amplitude, T_i is the period, ϕ_i is the phase, and m is the number of terms. The prior constraints on each sinusoidal term were assumed based on the periodogram of the white-light curve, while c₀ is assumed to follow $\mathcal{U}\left(0.5,1.5\right)$.

3.3.3 Reducing Common-Mode Systematics

A part of the systematics in spectroscopic light curves is instrument induced and does not vary with wavelength. This is the so-called “common-mode” systematics. Reducing the common modes when fitting spectroscopic light curves can lower the noise level and thus improve the S/N. One example is the divide-white method outlined in Gibson et al. (2013b) and Stevenson et al. (2014), where the common-mode systematics are approximated with the white-light systematics W(t) = f(t)/m*(t). The spectroscopic light curves are then divided by the white-light systematics, $$f_{\lambda }^{\ast }\left(t\right)=\frac{f_{\lambda }\left(t\right)}{w\left(t\right)},$$(11)

where f_λ(t) is the normalized spectroscopic light curve observed at wavelength λ. Since the broadband systematics can be be viewed as the average of narrowband systematics weighted by flux, W(t) is a good approximation to the common-mode systematics if the systematics in different narrowbands have weak wavelength dependences.

However, the spectroscopic systematics may be largely dependent on wavelength (or flux) if they are mainly induced by the astrophysical processes such as stellar activities or pulsations of either the target or the reference star. In this scenario, the divide-white method might incorrectly estimate the noise level in different wavebands, causing unexpected modification in some narrowband light curves. Therefore, in addition to the divide-white method, we also tested whether it would present a better performance if a scaling factor η_λ were added as a free parameter to account for the wavelength dependence of spectroscopic systematics. A similar method has been outlined in Mandell et al. (2013). We call this variant method as the divide-rescaled method, where the spectroscopic light curve is corrected to be $$f_{\lambda }^{\ast }\left(t\right)=\frac{f_{\lambda }\left(t\right)}{1+{\eta }_{\lambda }\left(w\left(t\right)-1\right)}.$$(12)

When η_λ = 1, the above equation is equivalent to Eq. (11), while η_λ approaches zero if the light curve at wavelength λ does not share the common-mode patterns.

In addition to these two methods, we also tried to fit the spectroscopic light curves directly without reducing the common-mode systematics, which is denoted as the direct-fit method. It serves as a null model such that we can investigate how a large amplitude of common-mode systematics would impact the transmission spectrum.

3.4 Bayesian Estimation and Model Comparison

The commonly used methods for Bayesian parameter estimation include the Markov chain Monte Carlo method (e.g., the Metropolis-Hastings algorithm, Metropolis et al. 1953; Hastings 1970; and the Hamiltonian Monte Carlo algorithm, Betancourt 2017) and nested sampling (Skilling 2004; Feroz & Hobson 2008). The MCMC method is highly effective for obtaining a stable estimate of the posterior joint distributions of the parameters. The nested sampling algorithm can also obtain the posterior estimates of the parameters, although it is initially designed for computing Bayesian evidence (also called the marginalized likelihood), which provides a criterion for model comparison. Since model comparison is a critical part in this work, we adopted the nested sampling method to obtain the Bayesian inferences. This method uses a set of sampling points to fully explore the prior space from the prior boundary toward higher-likelihood regions and is able to reveal multimodal distributions when implemented with the MULTINEST algorithm (Feroz et al. 2009). The sampling naturally comes to a convergence when the remaining parameter space has negligible evidence contributions. When the sampling is complete, we obtain both Bayesian evidence and the posterior estimates.

We used PyMultiNest (Buchner et al. 2014), which is built on the MULTINEST library, to estimate Bayesian evidence and parameter posteriors. In the context of Bayes’ rule, the model evidence $Z$ can be calculated by integrating the likelihood over the prior space of parameters Θ, which intrinsically is the average likelihood of a model hypothesis over its prior space. Following Feroz & Hobson (2008), this can be written as $$Z\left(\mathcal{D}|\mathscr{H}\right)={\displaystyle \int ℒ\left(\mathcal{D}|\Theta ,\mathscr{H}\right)}\pi \left(\Theta |\mathscr{H}\right){\text{d}}^D\Theta ,$$(13)

where Ɗ is the data, ℋ is the model hypothesis, ℒ is the likelihood function, π is the prior function, and D is the dimensionality of the parameter space. When the data are given, the model hypothesis with a lower dimensionality but a higher maximum likelihood usually presents stronger Bayesian evidence and thus prevails in model comparison.

In practice, we acquired ln $Z$ by calculating the natural logarithmic likelihood, $$\text{ln}ℒ\text{=}-\frac{N}{2}\text{lN}\left(2\pi \right)-\frac{1}{2}\text{ln}\left|K\right|-\frac{1}{2}{r}^T{K}^{-1}r,$$(14)

where N is the number of data points. The residual vector r and the covariance matrix K have been defined in Sect. 3.3.1 for the GP models. When using the BF models, the residual is r(t) = f(t) − m*(t)b(t), and K is just a diagonal matrix composed of jitter terms. We estimated ln $Z$ with 1000 sampling points, a sampling efficiency of 0.3, and an evidence tolerance of 0.1 when running PyMultiNest. The corresponding precision of ln $Z$ is 0.1–0.5, depending on model complexity.

We compared the model by comparing ln $Z$ between two models, which is equivalent to the natural logarithmic Bayes factor, $$\ln \frac{\text{pr}\left(\mathcal{D}|{\mathscr{H}}_1\right)}{\text{pr}\left(\mathcal{D}|{\mathscr{H}}_2\right)}=\ln \left(\frac{z_1}{z_2}\right),$$(15)

where the subscripts 1 and 2 denote hypothesis 1 (ℋ₁) and hypothesis 2(ℋ₂), respectively. We adopted the categories proposed by Kass & Raftery (1995) that ℋ₁ has very strong evidence against ℋ₂ if $\ln \left(Z_1/Z_2\right)\ge 5$;ℋ₁ has strong evidence against ℋ₂ if $3\le \ln \left(Z_1/Z_2\right)<5$; ℋ₁ has positive evidence against ℋ₂ if $1\le \ln \left(Z_1/Z_2\right)<3$; and two hypotheses have comparable evidence if | $\ln \left(Z_1/Z_2\right)<1$.

4 Results and Discussion

4.1 White-Light Curve Fitting

Figures 3 and 4 show the best-fit white-light curves of KIC 10657664 calibrated with four different reference stars, where the systematics are modeled with GPs and BFs, respectively. It is evident that the amplitudes of light-curve systematics vary with the adopted reference star. The light curve calibrated with REF-A1 exhibits the weakest systematics among the four groups of results, and its root-mean-squared error (RMSE) reaches ~1.3 times photon noise when fit with GP models. The other three light curves show stronger systematics with similar patterns. We speculate that the difference of systematics in amplitudes is mainly related to the distance between the star position and the rotation center of the OSIRIS field of view. As shown in Fig. 1, the three reference stars (REF-A2, REF-A3, and REF-A4) placed on CCD2 are farther away from the rotation center compared with the target star and the REF-A1 star. A similar scenario was reported in Nortmann et al. (2016), where the bump features in systematic noise might be attributed to vignetting in pupil space due to instrument rotation. Additionally, it is impossible to keep the five objects perfectly aligned with the central line of the long slit (7′ × 12″) during the observation. Therefore, discrepant flux variations might result from varying pixels or slight differential light losses produced by the finite slit width. We also tested whether the larger-amplitude systematics for REF-A2 to REF-A4 could be induced by differential sky variability, but found no correlation between them. In addition, the variability of the reference stars cannot explain the similarity of systematics in the light curves calibrated with REF-A2, REF-A3, and REF-A4 either.

The estimated transit parameters of KIC 10657664B are listed in Table A.2. We tested light-curve fitting assuming uninformative priors, but found that the transit parameters, particularly R₂/R₁, a/R₁, and i, failed to be constrained by the transit light curves using REF-A2, REF-A3, and REF-A4. Therefore we adopted Gaussian priors based on the estimates of Wong et al. (2020), which are well constrained by the collection of all 18 quarters of the Kepler long-cadence light curves. We used the transit depth obtained in the Kepler band (420–900 nm) to constrain that obtained in the OSIRIS R band (515–915nm) assuming that the transit depths of white dwarfs are independent of wavelength. Furthermore, our white-light curve analyses serve to extract precise and reliable broadband systematic noise, instead of improving the transit parameter estimation, which is critical for removing the common-mode systematics in the spectroscopic light curves.

The case of KIC 9164561 is more complicated. We used Gaussian priors to constrain R₂/R₁, a/R₁, and i based on the estimates of Rappaport et al. (2015) (18 Kepler quarters). Otherwise, the transit signals could not be correctly distinguished from the correlated noise through GP modeling because of the strong stellar pulsations and the short observation time out of transit. As described in Sect. 3.1, we additionally considered the effects of ELV and ILL when fitting the light curves of KIC 9164561. The best-fit white-light curves are shown in Fig. 5, and the corresponding transit parameters are listed in Table A.3. According to the results, we find that the M32 and the SHO kernels predict different systematics and central transit time ($t_c=-{0.00331}_{-0.00512}^{+0.00438}$ d when using the M32 kernel; $t_c=-{0.01426}_{-0.00157}^{+0.00171}$ d when using the SHO kernel; Table A.3). These discrepancies mainly arise because the boundaries of transit ingress and egress become indistinct due to the large amplitudes of oscillatory noise. We did not assume a Gaussian prior on t_c because the companion white dwarf exhibits a long-term transit timing variation reported by Rappaport et al. (2015), which might arise from an undiscovered third body with an orbital period of 8 to 14 years. Because the GTC transit observation is ~790 days (~623 binary orbits) later than the last transit observed by Kepler, it is difficult to tell which is closer to the true value, the M32-derived t_c or the SHO-derived t_c. Therefore, we kept both results derived from the two GP models and separately used their systematics for common-mode corrections in the subsequent spectroscopic light-curve analyses of KIC 9164561.

In addition to using GP models, we also tried to use various BF models to account for light-curve systematics. We conducted a Bayesian model comparison among 4096 different BFs for KIC 10657664 and eight sinusoidal models for KIC 9164561 based on their white-light curves. In the subsequent spectroscopic light curve fitting, we adopted those BFs with the strongest Bayesian evidence. The results of Bayesian evidence estimated by the nested sampling algorithm are listed in Table A.4 for KIC 10657664 and Table A.5 for KIC 9164561, where the evidence for GP models is also listed for comparison. Table A.4 shows that the best BF model changes with the adopted reference star, while all of the best models show correlations with the FWHM of the stellar PSF and the rotation angle of the instrument, suggesting the sources of light curve systematics. However, as shown in Fig. 4, the best BFs selected from a large number of models still have a poor performance in characterizing the light-curve systematics of REF-A3 and REF-A4. On the other hand, Table A.5 shows the best BF for KIC 9164561 consists of eight sinusoidal terms, but the corresponding fitting result is still worse than those of the GP models (Fig. 5). We note that the Bayesian evidence of the sinusoidal BF continues to increase with the number of harmonic terms, but we did not use more complex models with nine or more harmonic terms because their computational cost is high. For instance, when nine harmonic terms are used in Eq. (10), it takes about 10 h for a 64-core parallel computation (3 GHz CPUs) to fit one light curve. Meanwhile, adding more trivial terms with higher frequencies will not significantly improve the results. Therefore, we adopted m = 8 in Eq. (10) in the subsequent light-curve analyses of KIC 9164561.

4.2 Comparison of Transmission Spectra

The raw spectroscopic light curves of KIC 10657664 and KIC 9164561 are shown in Fig. A.2. We used the methods described in Sect. 3 to compute the transmission spectra of the two targets. We assumed uniform priors of U(0.03,0.3) on narrowband transit depths when fitting the spectroscopic light curves. We did not assume Gaussian priors on R₂/R₁ in spectroscopic light-curve fitting so as to avoid the resulting transmission spectra being largely affected by the informative priors. We compared the computed transmission spectra with the broadband transit depths to check the general offsets of transit depth. We analyzed the potential wavelength dependency of a transmission spectrum by fitting the spectrum with polynomials from the zeroth order up to the 19th order and then selected the best-fit polynomial with the strongest model evidence. A transmission spectrum is determined to be flat and featureless if the zeroth-order polynomial (i.e., a constant function) presents the strongest evidence. Otherwise, the transmission spectrum is determined to exhibit spurious spectral features resulting from data analysis procedures or noise fluctuation.

Fig. 3 Results of white-light curve fitting for KIC 10657664 using the GP-M32 model. The four split parts correspond to the light curves calibrated with four reference stars (REF-A1 to A4). In each part, the top panel shows the white-light curve (black dots), the best-fit model (red), and the GP systematics (green), the middle panel shows the detrended light curve (black dots) and the transit model (red), and the bottom panel shows the residuals, where εph denotes the uncertainty of photon-dominated noise. The green shaded intervals indicate 5σ confidence of the GP models. The GP systematics are shifted to match the light-curve baselines.

4.2.1 KIC 10657664B

Figure 6 shows the transmission spectra of KIC 10657664B computed with different reference stars (REF-A1 to A4), different noise models (the GP-M32 models and the best-fit BF models listed in Table A.4), and different common-mode removal techniques (direct-fit, divide-white, and divide-rescaled). Although most of the derived spectra can be considered to be flat and featureless, only the result corresponding to REF-A1, GP, and direct-fit has arelatively good precision and is consistent with the broadband transit depth. The statistical metrics for the spectra in Fig. 6 are listed in Table A.6.

Different reference stars. The transmission spectra corresponding to REF-A1, whose light curves have weaker systematics, are closer to the broadband transit depth and have stronger evidence to be featureless than the results from the remaining reference stars. This indicates that the amplitudes of systematics in spectroscopic light curves are critical for obtaining an accurate transmission spectrum. Using different reference stars may cause considerable discrepancies in transmission spectra due to different levels and trends of light-curve systematics. We note that it is commonly impractical to simultaneously monitor multiple reference stars when using long-slit spectroscopy. A reference star should be selected before the observation based on its physical properties, including the stellar type, magnitude, and light variation. Thus by comparing the results from different reference stars, we do not strive to select a proper reference star based on the resulting transmission spectrum. Instead, we use the different reference stars to represent the inconsistent instrument-specific systematics that would exist in different observations. We show that such discrepancies in light-curve systematics can significantly impact the results.

Different noise models. In Fig. 6, the GP models and the BF models lead to totally different results, indicating that the adopted noise model also has a large impact on the transmission spectrum. When using the GP model to fit the light-curve systematics, there is no significant transit depth offset in the derived transmission spectra. When the BF models are used, however, the transit depth offsets are quite large when the light-curve systematics are too strong to be corrected by the BF model (see the results using REF-A3 and REF-A4 in Fig. 6). This is mainly because in spectroscopic light-curve fitting, the orbital parameters (semimajor axis, inclination, mid-transit epoch) are fixed to the best-fit values from broadband light-curve fitting, which gives critical constraints on the transit model. However, if the adopted best-fit orbital parameters are biased from true values, the derived spectra may present considerable offsets. When GP models are used to fit the white-light curves, the derived orbital parameters, especially the mid-transit epoch, are basically consistent among four groups of data. When the BF models are used, however, the best-fit orbital parameters are less consistent. Therefore, the derived spectra using the GP models are closer to the correct transit depth than those using the BF models, as shown in Fig. 6. Another reason for the results of BF models being less stable is that we adopted the best BFs based on the white-light curve analysis and used them to fit spectroscopic systematics. Since the deterministic BF models are less flexible than stochastic GP models, it is difficult for the best BF model to consistently solve the systematic noise in all narrow-bands, and thus result in biased estimation on chromatic transit depths.

Different common-mode removal methods. According to Fig. 6, the common-mode removal methods in spectroscopic light-curve fitting also have large impacts on transmission spectra. When GP models are used, the spectra computed from the direct-fit method can have much larger errors than those from the divide-white and divide-rescaled methods. We note that the GP model usually captures the smooth variations in systematic noise. Therefore, when the common modes are not reduced, the GP model will attribute the jitter patterns induced by seeing variations to uncorrelated noise (i.e., white noise), which significantly reduces the light-curve S/N and thus overestimates the transit depth errors. Such large uncertainties will negate any possible detection from an exoplanet atmosphere. When the common-modes are reduced, however, the uncertainties of transmission spectra become so small that the results may show spurious features. When BF models are used, the common modes and the remaining components in light-curve systematics are fit by the same BF, resulting in similar uncertainties. As previously mentioned, the general offset of a transmission spectrum can be extremely large when the BF model fails to characterize the light-curve systematics. We find that reducing the common-mode patterns can slightly reduce such an offset, but does not solve the problem. This is because the best BF selected for the broadband systematics is too complex for the remaining spectroscopic systematics after the common-mode removal, which is likely to cause overfitted models and weaken the transit signals. In our preliminary test, when using the polynomial BF up to the first order in spectroscopic light-curve fitting, there is no significant general offset in the derived transmission spectrum, but in some wavebands, the light curve cannot be correctly fit. An ideal solution is to conduct a separate BF model selection for each narrowband light curve. We did not consider this method, however, because its computational cost is high. When we compare the divide-white and the divide-rescaled methods, their results are similar. The two methods have a risk of underestimating the transit depth errors and thus result in spurious features in the derived transmission spectra. Therefore it may be unnecessary to remove the common-mode noise when there are only small amplitudes of systematics in the light curves (e.g., when fitting the light curves calibrated with REF-A1).

Fig. 4 Same as Fig. 3, but the systematics are modeled with the best-fit BFs (#1) listed in Table A.4.

Fig. 5 Results of white-light curve fitting for KIC 9164561 using different noise models. Left column: GP-M32 model. Middle column: GP-SHO model. Right column: sinusoidal BF model. Top row: white-light curve (black dots), best-fit model (red), and systematic noise (green). Middle row: detrended light curves (black dots) and transit model (red). Bottom row: residuals. The green shaded intervals indicate the 5cr confidence of the GP models. The GP systematics are shifted to match the light-curve baselines. The systematics exhibit positive offsets due to the ILL and ELV effects.

Fig. 6 Transmission spectra of KIC 10657664B calculated using different noise models (left column: GP-M32 model. Right column: best-fit BF model), different reference stars (REF-A1 to A4 from the top row to the bottom row), and different common-mode removal techniques (direct-fit: blue circles; divide-white: green squares; divide-rescaled: red diamonds). The solid lines and shaded areas are the best-fit polynomials with the strongest Bayesian evidence and their 1er credible intervals, respectively. The black horizontal lines are the broadband radius ratio from the Kepler data. The data points are slightly shifted along the wavelength for clarity.

4.2.2 KIC 9164561Β

In Fig. 7, we compare the transmission spectra of KIC 9164561B separately computed with two GP models (GP-M32 and GP-SHO) and one sinusoidal BF model (m = 8), whose statistical metrics are listed in Table A.6. Similar to the previous target, there are considerable discrepancies in the spectra derived from different noise models, among which the spectrum calculated with the GP-SHO model and the divide-rescaled method is better than other spectra considering its median uncertainty and general offsets. When the M32 kernel is used, the resulting transmission spectra have larger uncertainties, especially when the common modes are not reduced. This suggests that the M32 kernel, along with other radial basis functions, has a poor performance in extracting quasi-periodic correlated noise. When the S HO kernel is used, the derived transmission spectra have a better precision and are consistent with the broadband transit depth. In contrast, we fail to recover correct transmission spectra when the sinusoidal models are used. We speculate that this is because the sinusoidal baseline model mainly characterizes the oscillatory noise induced by stellar pulsation, but barely includes the additional patterns induced by light-curve calibration, seeing variation, and other instrumental systematics. Therefore the estimated transit depth can be significantly biased in some wavebands. Reducing the common-mode noise helps to reduce these biases, but the consequent underestimation of transit depth errors still causes spurious features in transmission spectra. The results from the divide-rescaled method are better than those from the divide-white method because the main component of systematic noise originates from stellar pulsation, which strongly depends on wavelength. The divide-white method would cause excess or deficient corrections on the pulsation-induced systematics and is likely to result in false linear trends when GP models are used.

The example of KIC 9164561B indicates that the derived transmission spectra can be highly model dependent when the light-curve systematics exhibit oscillatory features and have an amplitude comparable with the transit signal. The oscillatory noise may blur the boundaries of transit ingress and egress and makes it impossible to constrain the limb-darkening features, thus impacting the estimated transit parameters. Even for flexible GP models, the results can be biased if the systematics are not characterized by suitable GPs. Therefore, it is indeed difficult to recover a transit light curve contaminated by strong stellar variations. Under these circumstances, we recommend using the SHO kernel as a robust model of stellar variations (Foreman-Mackey et al. 2017) before we have a better understanding on the noise sources and proper modeling on their variability.

4.3 Further Discussion of the Flux Dilution Correction

In Sect. 4.3, we consider the flux dilution effect in transit light curves due to the strong thermal emission of hot white dwarfs. We correct for the flux dilution by calculating the flux ratios of the binaries assuming pure blackbody radiation because there is no secondary eclipse observation to directly constrain the spectroscopic flux ratio. When we use fixed effective temperatures in Eq. (2), the derived flux ratio is proportional to the transit depth. However, the relation between flux ratio and radius ratio would be broken if the true spectroscopic flux ratio were to deviate from the blackbody hypothesis. Therefore, here we briefly discuss whether the biases of this correction would significantly impact the transmission spectra.

The absolute deviation of radius ratio δ_ℛ resulting from that of flux ratio ${\delta }_F$ is $${\delta }_{ℛ}=\frac{{ℛ}_{\lambda }}{2\left(1+{ℱ}_{\lambda }\right)}\delta ℱ,$$(16)

where ℛ_λ and ${ℱ}_{\lambda }$ denote the true values of radius ratio R₂/R₁ and flux ratio F₂/F₁ at wavelength λ. A brief derivation of Eq. (16) is presented in Appendix B.

Since the companion-to-host flux ratios for the binaries of KIC 10657664 and KIC 9164561 are measured to be ~0.0189 and ~0.0285 (Wong et al. 2020; Rappaport et al. 2015) in the Kepler band and are far lower than 1, the relative deviation of radius ratio δ_ℛ/ℛ_λ then approximately equals ${\delta }_{ℱ}/2$ according to Eq. (16). Thus even if we did not correct the flux dilution effect and made a huge error on flux ratio that ${\delta }_{ℱ}={ℱ}_{\lambda }\approx 0.02$, the corresponding relative deviation of radius ratio would be merely δ_ℛ/ℛ_λ The transmission spectra in Figs. 6 and 7 have relative errors of approximately 2.5% for KIC 10657664B and 1.8% for KIC 9164561B when they are calculated with the GP-M32 model and the divide-white method. The flux dilution correction can only introduce very limited deviations that are smaller than the 1er credible intervals in the radius ratio. Therefore, we conclude that the significant spectral features or general offsets in the derived transmission spectra should not arise from the flux dilution effect.

Fig. 7 Transmission spectra of KIC 9164561B calculated with different noise models (Top panel: GP-M32 model. Middle panel: GP-SHO model. Bottom panel: sinusoidal baseline model) and different common-mode removal methods (direct-fit: blue circles; divide-white: green squares; divide-rescaled: red diamonds). The solid lines and shaded areas are the best-fit polynomial trends with the strongest Bayesian evidence and their 1σ credible intervals, respectively. The black lines are the broadband radius ratio from the Kepler data. The data points are slightly shifted along the wavelength for clarity.

5 Conclusions

We used two transiting white dwarfs, KIC 10657664B and KIC 9164561B, as hot-Jupiter analogs to test the methods that are commonly used in low-resolution transmission spectroscopy. This benchmark test allowed us to compare the outcomes of different data analysis techniques. Its advantage is that we know beforehand that a white dwarf transmission spectrum should be flat and featureless. The transit events of these two targets were observed with the GTC OSIRIS instrument. The acquired data were reduced and analyzed in the same way as exoplanet transit spectroscopic observations.

In the analyses of the two targets, we obtained discrepant transmission spectra when we used different noise models. The results show that GP models have a better robustness than BF models when there are large amplitudes of systematic noise. However, different GP kernels may result in different transit parameters and final transmission spectra, especially when there are strong oscillatory noises in the light curves. Thus, the most suitable model must be carefully chosen based on model comparison.

In the example of KIC 10657664, the different reference stars result in different intensities of light-curve systematics that rely on specific instruments and star positions. These differences can cause significant discrepancies in the derived transmission spectra. In the particular case of GTC OSIRIS, we suggest that it is better to select a reference star closer to the target star if there are several options, given consideration of their physical properties including the spectral type, brightness, and light variations. In the example of KIC 9164561B, the GP-SHO model is found to have stronger Bayesian evidence and produces more accurate transmission spectra than the GP-M32 and the sinusoidal BF models in characterizing the stellar pulsation noise. We expect, however, that similar situations will apply to other ground-based observations of exoplanets, and the resulting large uncertainties could mute any weak exoplanet atmospheric signals.

We showed that reducing the common-mode systematics with the divide-white or the divide-rescaled method helps to reduce the spectroscopic systematic noise, but also brings in the risk of underestimating the transit depth errors and introducing additional spurious spectral signatures. Thus it may be unnecessary to conduct a common-mode removal if there is only a small amplitude of systematics.

In summary, we conclude that no model exists that would always ensure a perfect characterization of the systematic noise for ground-based low-resolution transit observations. The bias of extracted systematics inevitably causes the bias of transmission spectra. Therefore, when transmission spectroscopy studies are conducted, it is necessary to validate the results based on various model comparison, using either simulated data or real data, such as those from transiting white dwarfs, to determine their model dependence. In particular, small wiggles within the transit depth errors in a transmission spectrum should be interpreted with caution. Furthermore, the target of interest should ideally be repeatedly observed with different instruments to reduce the impact of instrumental systematics and increase the reliability of the results.

Acknowledgements

The authors thank the anonymous referee for their valuable comments and suggestions. G.C. acknowledges the support by the B-type Strategic Priority Program of the Chinese Academy of Sciences (Grant No. XDB41000000), the National Natural Science Foundation of China (Grant Nos. 42075122, 12122308), the Natural Science Foundation of Jiangsu Province (Grant No. BK20190110), Youth Innovation Promotion Association CAS (2021315), the China Manned Space Project (CMS-CSST-2021-B12), and the Minor Planet Foundation of the Purple Mountain Observatory. This work is based on the observations made with the Gran Telescopio Canarias installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma.

Appendix A Additional Tables and Figures

Fig. A.1 Time series of the state vectors used in the BF models for KIC 10657664. From top to bottom: Target position drift in the dispersion direction, target position drift in the spatial direction, the FWHM of the target PSF in the spatial direction, the airmass, and the instrumental rotation angle.

Fig. A.2 Raw spectroscopic light curves of KIC 10657664 and KIC 9164561. The colors of the light curves indicate corresponding wavebands. The light curves are arbitrarily shifted for clarity. The flux errors are too small to be shown.

Appendix B: Derivations of Eq. 16

We denote the true values of the companion-to-host radius ratio and flux ratio as ℛ and $ℱ$, respectively. The flux dilution correction in Eq. 1 is equivalent to $$1-m^{\ast }\left(t\right)=\frac{1-m\left(t\right)}{1+ℱ},$$(B.1)

where m*(t) is the diluted transit model, while m(t) is the model without the dilution effect.

At the mid-transit epoch t_c, there is a linear approximation between the normalized flux and the transit depth $$\begin{array}{lll}m{|}_{t=t_c}\hfill & =\hfill & 1-\xi {ℛ}^2{+}_o\left({ℛ}^2\right),\hfill \\ m\ast {|}_{t=t_c}\hfill & =\hfill & 1-\xi {ℛ}^{\ast 2}{+}_o\left({ℛ}^{\ast 2}\right),\hfill \end{array}$$(B.2)

where ξ is a constant accounting for the host star’s limbdarkening effect, and ℛ* is the diluted radius ratio. Substituting Eqs. B.2 into Eq. B.1, we eliminate ξ and obtain $$ℛ=\sqrt{1+ℱ{ℛ}^{\ast }}.$$(B.3)

When error propagation is considered, the deviation of radius ratio δ_ℛ arising from that of flux ratio ${\delta }_{ℱ}$ writes $$\frac{{\delta }_{ℛ}}{ℛ}\approx 0.5{\delta }_{ℱ}=\frac{ℛ}{2\left(1+ℱ\right)}{\delta }_{{}_{ℱ}},$$(B.4)

where ℛ* is eliminated. Thus for $0<ℱ\ll 1$, the relative deviation of radius ratio is $$\frac{{\delta }_{ℛ}}{ℛ}\approx 0.5{\delta }_{ℱ}.$$(B.5)

References

All Tables

All Figures

	Fig. 1 Through-slit images of the two observations. Upper row: Positions of KIC 10657664 and the reference stars (REF-A1 to A4) in CCD1 and CCD2. Lower row: Positions of KIC 9164561 and the reference star (REF-B1) in CCD1. The red mark is the GTC Nasmyth rotator center. Only a part of the data section is displayed for clarity.
In the text

	Fig. 2 Time-averaged stellar spectra of all objects. The upper panel shows the spectra of KIC 10657664 and its four reference stars. The lower panel shows the spectra of KIC 9164561 and its single reference star. The green shaded intervals indicate the wavebands for spectroscopic light curves.
In the text

Fig. 3 Results of white-light curve fitting for KIC 10657664 using the GP-M32 model. The four split parts correspond to the light curves calibrated with four reference stars (REF-A1 to A4). In each part, the top panel shows the white-light curve (black dots), the best-fit model (red), and the GP systematics (green), the middle panel shows the detrended light curve (black dots) and the transit model (red), and the bottom panel shows the residuals, where εph denotes the uncertainty of photon-dominated noise. The green shaded intervals indicate 5σ confidence of the GP models. The GP systematics are shifted to match the light-curve baselines.

In the text

	Fig. 4 Same as Fig. 3, but the systematics are modeled with the best-fit BFs (#1) listed in Table A.4.
In the text

Fig. 5 Results of white-light curve fitting for KIC 9164561 using different noise models. Left column: GP-M32 model. Middle column: GP-SHO model. Right column: sinusoidal BF model. Top row: white-light curve (black dots), best-fit model (red), and systematic noise (green). Middle row: detrended light curves (black dots) and transit model (red). Bottom row: residuals. The green shaded intervals indicate the 5cr confidence of the GP models. The GP systematics are shifted to match the light-curve baselines. The systematics exhibit positive offsets due to the ILL and ELV effects.

In the text

Fig. 6 Transmission spectra of KIC 10657664B calculated using different noise models (left column: GP-M32 model. Right column: best-fit BF model), different reference stars (REF-A1 to A4 from the top row to the bottom row), and different common-mode removal techniques (direct-fit: blue circles; divide-white: green squares; divide-rescaled: red diamonds). The solid lines and shaded areas are the best-fit polynomials with the strongest Bayesian evidence and their 1er credible intervals, respectively. The black horizontal lines are the broadband radius ratio from the Kepler data. The data points are slightly shifted along the wavelength for clarity.

In the text

Fig. 7 Transmission spectra of KIC 9164561B calculated with different noise models (Top panel: GP-M32 model. Middle panel: GP-SHO model. Bottom panel: sinusoidal baseline model) and different common-mode removal methods (direct-fit: blue circles; divide-white: green squares; divide-rescaled: red diamonds). The solid lines and shaded areas are the best-fit polynomial trends with the strongest Bayesian evidence and their 1σ credible intervals, respectively. The black lines are the broadband radius ratio from the Kepler data. The data points are slightly shifted along the wavelength for clarity.

In the text

	Fig. A.1 Time series of the state vectors used in the BF models for KIC 10657664. From top to bottom: Target position drift in the dispersion direction, target position drift in the spatial direction, the FWHM of the target PSF in the spatial direction, the airmass, and the instrumental rotation angle.
In the text

	Fig. A.2 Raw spectroscopic light curves of KIC 10657664 and KIC 9164561. The colors of the light curves indicate corresponding wavebands. The light curves are arbitrarily shifted for clarity. The flux errors are too small to be shown.
In the text