Constraining the oblateness of transiting planets with photometry and spectroscopy

Previous studies compared the transit light curve of an oblate planet to that of a spherical planet and showed that the oblateness signal manifests itself at the ingress and egress phases (e.g., Carter & Winn 2010a; Zhu et al. 2014). The oblateness-induced signal is a geometrical effect and is obtained as the residuals from fitting the transit observation of an oblate planet with a spherical planet model. The signal arises mostly due to difference in contact times at the stellar limb (ingress and egress) between the transiting oblate planet and the corresponding spherical planet of the same area. Since the geometry of a transiting planet is the same when taking photometric and spectroscopic transit measurements, the oblateness signature should also manifest itself in the transit RM signal.

To illustrate the oblateness signature in the light curve and RM signal, we follow the work of Carter & Winn (2010a) which analysed seven Spitzer transit observations of the giant planet HD 189733b to constrain its oblateness. The planet was selected due to its large size, bright host star and availability of precise Spitzer data. First, they calculated the theoretical oblateness signal amplitude expected for HD 189733b by simulating its photometric light curve assuming Saturn-like oblateness of $f=0.098$ and projected obliquity $\theta=45\degree$. The planet has the following parameters: planet-star mean radius ratio $\bar{R}_{p}=0.15463\,R_{\ast}$ (equivalent to $1.13\,R_{\mathrm{Jup}}$), scaled semi-major axis $a/R_{\ast}=8.81$, inclination $i_{p}=85.58\degree$ (or impact parameter $b=0.68$) and orbital period $P_{\mathrm{orb}}=2.22$ days (Torres et al., 2008). The host star is of K2 spectral type with V-magnitude $m_{V}$ = 7.7, radius $R_{\ast}=0.75\,R_{\odot}$ and $v\sin{i_{\ast}}=3.5$ km s^-1 but we assume a stellar inclination $i_{\ast}=90\degree$ (rotation axis parallel to sky plane). We also simulated mock oblate planet transit signals (light curve and RM signal) for this planet using SOAP3.0. The light curve was generated with Spitzer cadence of 2 minutes using quadratic limb darkening coefficients (LDCs) of [0.076, 0.034] in the Spitzer infrared band as given in Carter & Winn (2010a). The RM signal was generated with an integration time of 4 mins²²2 This integration time is similar to the 5 mins exposures of HARPS archival RM observations of HD 189733b used by Triaud et al. 2009; Wyttenbach et al. 2015. which, for real observations, will help reduce the amplitude of stellar noise in this star as recommended in Chaplin et al. (2019) for K-dwarfs. We obtained LDCs for the star in the visible band using Claret & Bloemen (2011) and assumed $\lambda=0\degree$ (close to the value of $0.85\degree$ derived in Triaud et al. 2009). We then separately fit the oblate planet light curve and also its RM signal with spherical planet models. The residuals (oblate - spherical) from both fits are shown in Figure 2 with their amplitudes ($S_{\mathrm{obl}}$) quoted. The residuals show the effects of oblateness at ingress and egress as described in Seager & Hui (2002) and Barnes & Fortney (2003).

The residuals from our light curve fit is in agreement, in shape and amplitude, with the photometric result obtained in Carter & Winn (2010a) for this planet (see their Figure 1). The oblateness-induced signature in the RM signal, shown here for the first time, has an amplitude $S_{\mathrm{obl}}=0.75$ m s^-1 whereas the ESPRESSO spectrograph (installed on the 8m class telescopes of the ESO VLT; Pepe et al. 2014) is capable of attaining instrumental RV precisions up to 0.1 m s^-1. Using the ESPRESSO Exposure Time Calculator (ETC)³³3www.eso.org/observing/etc - Version P105.6, we estimated that the theoretical ESPRESSO RV precision for this star for the 4-minute exposure time is 0.31 m s^-1 for observation in the high resolution (UT1; 140,000) mode⁴⁴4Note: In simulating the aforementioned RM signal, we used this ESPRESSO resolution to calculate the FWHM of the CCF which was additionally broadened to account for macro-turbulence (convection) using calibration equation from Doyle et al. (2014).. Taking into account the rotational velocity broadening which degrades the quality factor of the spectra with respect to its non-rotating counterpart, the RV precision for this star degrades by a factor of 1.8 to $\sim$0.56 m s^-1 (calculated by interpolating between values in Table 2 of Bouchy et al. 2001). This precision can still be better if data is phase folded over an increased number of transit observations (n) as it scales in the photon noise limit with $1/\sqrt{n}$. If we assume that seven transits are observed (number of Spitzer transits analysed in Carter & Winn 2010a), a better RV precision of $\sim$0.2 m s^-1 would be achieved which is much lower than the spectroscopic oblateness amplitude of this system. Moreover, RM measurements of longer period planets can be obtained with longer exposure times that provide a higher RV precision such that several transit observations are not necessarily needed.

Although the analysis of this system is simply illustrative, it hints that ESPRESSO RM measurements of "relevant" systems with transiting planets could allow reasonable constraints to be placed on planet oblateness in addition to those from light curve analysis. We further investigate this possibility in the following sections with suitable considerations.

We point out that the oblateness signal can be similar to the signal caused by the presence of rings (Seager & Hui, 2002; Barnes & Fortney, 2004). However, depending on ring orientation and size of the gap between the planet and ring, the ring signature can be distinct showing more peaks at ingress and egress than the oblateness signature (see Ohta et al. 2009; Akinsanmi et al. 2018). Additionally, an estimate of the planetary density can help distinguish them since a planet with rings will appear larger and cause the planetary density to be underestimated (Zuluaga et al., 2015; Akinsanmi et al., 2020).

In the previous section we showed the signature of oblateness in the photometric light curve and spectroscopic RM signal for a short period planet. However, short period planets such as HD 189733b have short tidal evolution timescales so they are expected to already have circularized orbits and synchronised rotations such that $P_{rot}$ becomes equal to $P_{\mathrm{orb}}$ (Guillot et al., 1996). Such short period planets cannot have significant rotation induced oblateness since their rotations will be too slow (Seager & Hui, 2002). Indeed, Carter & Winn (2010a) did not detect oblateness in this planet for this same reason inspite of the large theoretically expected photometric signal (top panel of Fig. 2). However, short period planets can still be radially deformed (in direction of the star) due to the strong stellar tidal interaction. Previous studies (e.g. Correia 2014; Akinsanmi et al. 2019; Hellard et al. 2019) discussed the detectability of tidally induced deformation in short period planets and showed that this deformation diminishes with distance to the star.

Strong tidal interaction also affects the obliquity of short period planets by driving the value of $\theta$ to zero at the same short timescale for attaining rotation synchronization (Peale, 1999). Avoiding stellar tidal interaction then imposes the requirement for long period orbits to ensure rapid planet rotation for significant oblateness. Seager & Hui (2002) and Carter & Winn (2010b) showed that a Jupiter-mass planet orbiting a Solar twin star at a distance of $\sim$0.2 AU ($P_{\mathrm{orb}}\simeq 30$ days) will have tidal synchronization timescale of $\sim$10 Gyrs (which is greater than the age of most host stars) so that the planet is not expected to have spun down and can therefore have significant rotation.

The spin axis of an oblate planet will precess due to the gravitational torque exerted on the planet by the host star. The effect of spin precession on transit signals is that the orientation of the oblate planet changes with time and so does its projection. This leads to transit depth variations between obtained light curves and amplitude variations between the RM signals. Such variations can complicate efforts to measure the ingress and egress oblateness signature since combining successive transits might average out the subtle signal. However, a Jupiter or Saturn-like planet with $P_{\mathrm{orb}}\simeq 30$ days is expected to have precession period of $\sim$50 yrs (Carter & Winn, 2010b) which is too long to be observed within a few transit observations of the planet. For example, the spin axis of a $P_{\mathrm{orb}}=30$ days planet will precess by a negligible $1.8\degree$ in 3 successive transit observations thereby allowing the phase folding of data to probe oblateness.

Although some studies (e.g. Carter & Winn 2010b; Biersteker & Schlichting 2017) have illustrated the possibility of detecting oblateness due to spin precession for $P_{\mathrm{orb}}<30$ days, we investigate here the case of planets with $P_{\mathrm{orb}}\geq 30$ days for which significant rotation induced oblateness is expected and spin precession is negligible.

As the oblateness signal can depend strongly on the projected orientation of the planet ($\theta$) and that of its orbit ($i_{p}$ or $b$), it is important to identify the combinations of these two parameters that are optimal for detecting oblateness. Even though the value of $\theta$ is unknown a-priori for planets, an understanding of how it affects the detectability of oblateness is crucial as it can be used to calculate the maximum theoretical oblateness signal expected for a planet and thus aid in target selection. Although previous works (e.g., Barnes & Fortney 2003; Zhu et al. 2014) have identified $\theta-b$ parameter combinations with the highest photometric oblateness signal amplitude, we are particularly interested in understanding how the spectroscopic signal amplitude varies in comparison to the precision of new spectrographs like ESPRESSO.

We take case study of the HD189733 system assuming its Jupiter-sized planet ($1.13\,R_{\mathrm{Jup}}$) has a circular orbit around the star with a longer period of 50 days (based on requirements from previous section). As in the previous section, we assume a projected planet oblateness of $f=0.098$. We also assume that the star rotates slightly faster with $v\sin{i_{\ast}}=5$ km s^-1 (taking $i_{\ast}=90\degree$). We use quadratic LDCs for this star in the visual band obtained from Claret & Bloemen (2011) but re-parameterized as [$q_{1}$, $q_{2}$] following Kipping (2013). The full adopted parameters for the simulated system is given in Table 1.

We used SOAP3.0 to generate theoretical transit signals (light curves and RM signals) for the adopted oblate planet on a grid of obliquity ($\theta$) and impact parameter ($b$) values with a total of 42 combinations (excluding grazing transits since their ingress/egress are undefined). The transit durations and ingress (egress) durations of the transit signals vary with $b$. From $b$ = 0 to $b$ = 0.8, the transit durations decrease from $\sim$6.3 hrs to $\sim$4.5 hrs whereas the ingress (egress) durations increase from $\sim$55 mins to $\sim$90 mins. Transits at higher impact parameters will thus allow better sampling of oblateness induced features due to their longer ingress and egress durations. The light curves were simulated with 2-min integration time (or binning) similar to TESS but applicable to other photometric instruments to increase the attained precision of each measurement. The RM signals had a longer integration time of 8 mins⁵⁵5 Spectroscopic transit observations usually require longer integrations than the photometric because spectrographs lose photons due to slit losses, stray light, and scattered light and so require longer time to collect more photons (Oshagh, 2018a). Long integrations are also used to average out and reduce stellar RV noise. We note also that longer period planets can have longer integrations which should be chosen such that the ingress and egress are still well sampled. which enables ESPRESSO reach a higher RV precision of 0.22 m s^-1 while still providing good temporal resolution of the ingress and egress of this planet. Other spectrographs (like HARPS installed on a smaller 3.6 m telescope) are not capable of reaching such precision within this short integration time. To attain the same precision as ESPRESSO, HARPS will require $\sim$30 mins integrations which will not allow enough data points to probe oblateness signatures at ingress and egress.

We perform a least squares fit to the oblate planet light curve and RM signal of each parameter combination using a spherical planet model and obtain the amplitude, $S_{\mathrm{obl}}$, of the residuals (at ingress/egress). The free parameters in the light curve fit were $R_{p}$, $a/R_{\ast}$, $b$, $q_{1}$, and, $q_{2}$ while the RM signal fit additionally had $v\sin{i_{\ast}}$ and $\lambda$.

Figure 3 shows the contour plots generated using the residual amplitudes from fitting the photometric light curve and spectroscopic RM signal of each parameter combination. It represents the amplitude of observable oblateness-induced signal, $S_{\mathrm{obl}}$, at each parameter combination. The contours are shown only for positive $\theta$ angles as the pattern is symmetric about $\theta=0$. In both contour plots, we see a similar trend showing that the amplitude of photometric and spectroscopic oblateness signal is lowest (yellow regions) at zero obliquity across all impact parameters. This implies that the best-fit spherical model can easily emulate the light curve and RM signal of the oblate planet by adjusting its fit parameters thereby making it difficult to detect oblateness at these orientations. However, at higher obliquities, the amplitude increases with impact parameter due to asymmetry between ingress and egress of the oblate planet transit signal that cannot be easily emulated by a spherical model. The oblateness signal reaches its peak amplitude at points around $\theta=45\degree$, $b=0.7$ (dark blue regions) in agreement with the results of earlier photometric studies (e.g. Seager & Hui 2002; Carter & Winn 2010a) but confirmed here to be same in spectroscopy. The contour plot is nearly symmetric about the vertical $45\degree$ line with the exception of a few orientations around $b,\,\theta=0,\,90\degree$.

We also investigated how the spin-orbit angle ($\lambda$) affects the amplitude of oblateness signal. For the yellow regions in Fig. 3 with low oblateness signal amplitudes (e.g. orientations with $\theta\approx 0\degree$ and those with $b\approx 0$), we find that the spectroscopic oblateness signal amplitudes can increase for planets with spin-orbit misalignment (i.e. $\lambda$ $\neq$ 0). For example, the spectroscopic oblateness amplitude at the orientation $\theta$ = $0\degree$, $b$ = 0.2 is only 0.05 m s^-1 when $\lambda=0\degree$ but significantly increases to 0.61 m s^-1 for $\lambda=30\degree$. In contrast, this orientation has a photometric oblateness amplitude of only 3 ppm irrespective of $\lambda$. The orientations in the yellow regions can thus be more favorable for detecting oblateness in spectroscopy than in photometry. Figure 8 shows the variation of the spectroscopic oblateness amplitude with $\lambda$ at orientations with $\theta=0\degree$ and also orientations with $b=0$.

In summary, Fig. 3 shows that the oblateness-induced signal is more prominent, in photometry and spectroscopy, for oblique planets at high impact parameters particularly for the points around $\theta=45\degree$, $b=0.7$ (optimal transit orientation) which provides the best opportunity to measure oblateness. The residual from the fit of this orientation is shown in Fig. 4 with photometric and spectroscopic oblateness amplitudes of $\sim$272 ppm and 1.1 m s^-1, respectively.

The amplitude of the photometric and spectroscopic oblateness-induced signal for a particular parameter combination ($\theta,\,b$) scales with the oblateness and planet radius as $\sim$ $(f/0.098)\,(\bar{R}_{p}/0.1546)^{2}$ (Barnes & Fortney, 2003). Therefore, planets with larger planet-to-star radius ratios will have larger oblateness signatures for the same projected oblateness. For the same planet, the spectroscopic oblateness amplitude ($S_{\mathrm{obl}}^{~{}~{}\mathrm{sp.}}$) additionally scales with the projected stellar rotational velocity following the relation

Detecting planet oblateness implies obtaining a measurement for the parameter $f$ by fitting a transit observation, of sufficient precision, with an oblate planet model. Recovering $f\simeq 0$ implies that a spherical planet model is a better fit to the observation than an oblate planet model. Detecting oblateness thus requires that $f$ is recovered with accuracy and statistical significance above zero from the fitting process.

To illustrate the photometric and spectroscopic detectability of oblateness, we simulated the transit signals (light curve and RM) of the oblate planet at the determined optimal orientation (Table 1 with $\theta\,=\,45\degree,\,b\,=\,0.7$) which have oblateness amplitudes of 272 ppm and 1.1 m s^-1 (Fig. 4). The light curve and RM signal were simulated with cadences of 2 mins and 8 mins, respectively. Random Gaussian (white) noise of different levels, $N$, was added to the simulated transit signals. We then investigated how well we can recover $f$ and at what noise level it would be difficult to distinguish between the oblate planet and spherical planet transit signals (light curve and RM). This is useful in order to understand the instrumental precisions required for photometric and spectroscopic detection of oblateness.

We performed a Markov Chain Monte Carlo (MCMC) analysis with SOAP3.0 model to recover the oblate planet transit parameters along with their uncertainties. We use the emcee package (Foreman-Mackey et al., 2013) with priors on the parameters as given in Table. 2. We assume that the LDCs can be known with 10% accuracy so we adopt gaussian priors on them centered on the true simulated values with 10% standard deviation. The MCMC was performed with 36 walkers each having 20,000 steps (which was several times the computed auto-correlation time as recommended in emcee as a convergence diagnostic). The initial 25% of the steps were then discarded as burn-in.

Figure 5 shows the oblateness detectability plot, in photometry and spectroscopy, indicating the median and standard deviations of the recovered $f$ at different noise levels (average of three MCMC realizations). For oblateness to be confidently detected ($f$ measured), we require that $f$ is recovered with $3\sigma$ significance above zero.⁶⁶6 As a check of our MCMC analysis, we confirmed that fitting a spherical planet transit signal with an oblate planet model recovers $f\simeq 0$ at different noise levels. We see in the spectroscopic detectability plot (right panel of Fig. 5) that $f$ is detected with $3\sigma$ significance for noise levels up to 1 m s^-1/8 mins. At higher noise levels, the distribution of the recovered $f$ samples include $f=0$ at $3\sigma$ implying that a spherical planet model is also probable. In the photometric detectability plot (left panel of Fig. 5), a $3\sigma$ detection is attained for noise levels up to 400 ppm/2 mins. We note that $2\sigma$ detection of oblateness can still be attained at higher noise levels (up to 550 ppm in photometry and 1.5 m s^-1 in spectroscopy).

Figure 6 shows realized posterior distributions of the retrieved parameters at the detectable noise limits (400 ppm noise added to light curve and 1 m s^-1 added to RM signal). In both cases, we see from the $f-\theta$ joint distributions that $\theta$ is not well-constrained and is strongly correlated with $f$ ($\theta$ is only better constrained for highly precise transit signals). Figure 3 already showed that the amplitude of the oblateness signal is fairly symmetric about $\theta=45\degree$ and reduces as the value of $\theta$ gets farther from $45\degree$. Therefore, for $\theta$ values different from $45\degree$ in the $f-\theta$ distribution, the amplitude of oblateness signal reduces and a higher value of $f$ is needed to fit the observation. The degeneracy between $f$ and $\theta$ is responsible for the long tail towards large oblateness in the $f$ distribution (also seen for all noise levels in Fig. 5).

We converted the posterior of the fitted mean radius $\bar{R}_{p}$ to $R_{eq}$ using equation 3 to show the evidently strong correlation between $f$ and $R_{eq}$; as the oblateness increases, the equatorial radius gets more elongated compared to the polar radius. It is interesting to see that the limb darkening parameters ($q_{1}$, $q_{2}$) are not strongly correlated with $f$ implying that very precise determination of the LDCs are not required to detect planetary oblateness which is contrary to the case for detecting rings (Akinsanmi et al., 2018). Indeed when we adopt a non-informative uniform prior of $\mathcal{U}(0,1)$ on the LDCs, $f$ is still similarly recovered but with larger uncertainties on the LDCs and on the other parameters. The posteriors from the RM signal MCMC (left panel in Fig. 6) also shows that $\lambda$ and $v\sin{i_{\ast}}$ are not correlated with $f$ implying that planetary oblateness does not affect the accurate retrieval of these parameters from RM signals.

We investigated also the retrieval of $f$ from combined analysis of the photometric light curve and spectroscopic RM signal. As seen in Fig. 5, $f$ is not recovered with $3\sigma$ for noise levels of 450 ppm (added to the light curve) and 1.2 m s^-1 (added to the RM signal). However, simultaneously fitting both transit signals allows the recovery of $f$ with $3\sigma$ significance and with better accuracy as shown in Fig. 7 (compared to the Figs. 5 and 6).