# On The Effects of Planetary Oblateness on Exoplanet Studies

David Berardo  Department of Physics and Kavli Institute for Astrophysics and
Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139,
USA  FRQNT Doctoral Research Scholarship  Julien de Wit  Department of Earth,
Atmospheric and Planetary Sciences, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA David Berardo<EMAIL_ADDRESS>
###### Abstract

When studying transiting exoplanets it is common to assume a spherical planet
shape. However short rotational periods can cause a planet to bulge at its
equator, as is the case with Saturn whose equatorial radius is almost 10%
larger than its polar radius. As a new generation of instruments comes online,
it is important to continually assess the underlying assumptions of models to
ensure robust and accurate inferences. We analyze bulk samples of known
transiting planets and calculate their expected signal strength if they were
to be oblate. We find that for noise levels below 100ppm, as many as 100
planets could have detectable oblateness. We also investigate the effects of
fitting spherical planet models to synthetic oblate lightcurves. We find that
this biases the retrieved parameters by several standard deviations for
oblateness values $>$ 0.1-0.2. When attempting to fit an oblateness model to
both spherical and oblate lightcurves, we find that the sensitivity of such
fits is correlated with both the SNR as well as the time sampling of the data,
which can mask the oblateness signal. For typical values of these quantities
for Kepler observations, it is difficult to rule out oblateness values less
than $\sim$0.25. This results in an accuracy wall of 10-15$\%$ for the density
of planets which may be oblate. Finally, we find that a precessing oblate
planet has the ability to mimic the signature of a long-period companion via
transit timing variations, inducing offsets at the level of 10s of seconds.

##  1 Introduction

Young planets which have yet to come into tidal equilibrium with their host
star may be rotating at speeds sufficient to cause their equators to bulge and
thus deviate substantially from a sphere. When analyzing the lightcurves of
transiting exoplanets it is typical to assume that a planet is perfectly
spherical, as is the case of the many fitting routines which implement the
transit models of Mandel & Agol (2002). Within our own solar system however,
we observe that the equatorial radii of Saturn and Jupiter are $9.8\%$ and
$6.5\%$ larger than their polar radii. Planets can become distorted due to
tidal forces from their host star or, as is the case with Jupiter and Saturn,
due to rapid rotation, causing them to bulge uniformly about their rotation
axis. The amount by which a planet is oblate is tied to its rotation rate as
well as its internal structure, which determines its deformability. The
ability to measure the oblateness of a planet would thus allow us to better
understand its current rotational properties and internal dynamics, as well as
its formation and evolution (which determine its current rotation) (Lissauer,
1995).

Figure 1: Left: A schematic showing the configuration of the planet and all
relevant quantities. The blue (polar) minor axis is the planet’s axis of
rotation, which is scaled down by a factor of $(1-f_{\mathrm{proj}})$ relative
to the red (equatorial) major axis. The planet may also be tilted by an angle
$\theta_{\mathrm{proj}}$ relative to its orbital plane. Note these quantities
are the sky-plane projected values an observer would see. Right: The
difference in the lightcurve of an oblate planet (solid) to a spherical planet
(dashed) with the same projected area. For a uniformly bright star, there is
no difference once the planets are fully inside the stellar disk. Note the
asymmetry induced by the orbit having a non-zero impact parameter.

The effect of oblateness on the lightcurve of a transiting planet has
previously described in the works of Seager & Hui (2002), Barnes & Fortney
(2003), Carter & Winn (2010a), and Carter & Winn (2010b). These works
highlight the features of an oblate lightcurve, which is primarily a variation
in the ingress and egress of their transits on the order of 10’s to 100’s of
ppm for the most optimal of planets. These properties outline the difficulties
in detecting oblateness through lightcurve variations, which has thus far
proven to be a difficult task. In Zhu et al. (2014), short cadence planets in
the Kepler survey were searched, yielding a tentative result of oblateness for
the planet Kepler-39b (KOI432.01). A further search was conducted in
Biersteker & Schlichting (2017), which attempted to detect oblateness using
depth variations across multiple transits as an oblate planet precesses over
time. This yielded yet another tentative result of oblateness for the warm
Saturn Kepler-427b. Recently Akinsanmi et al. (2020) studied the potential of
spectroscopic measurements to detect oblateness, finding that a combination of
photometry and spectroscopy has the potential to provide more accurate and
precise oblateness measurements.

In section 2 we outline the relevant quantities and processes that control
planetary oblateness. We then analyze a bulk sample of planets to determine
how many should have detectable oblateness for varying levels of sensitivity.
In section 3 we fit spherical planet models to synthetic oblate lightcurves to
study biases on the retrieved transit parameters. As shown in de Wit et al.
(2012), orbital parameters can compensate for shape and brightness
distribution effects on the shape of transit and eclipses, an effect which we
seek to further understand. In section 4, we fit oblate models to synthetic
and real planets in order to determine the regimes in which such models can
properly measure or rule out oblateness, as well as to determine if
populations can be separated into oblate/non-oblate based on properties such
as orbital period. Finally in section 5 we discuss the effects oblateness may
have on the precision of measured planet densities, as well as the ability for
time-varying oblateness to mimic transit timing variations due to precession
of a planet’s orbital axis.

##  2 Description and Physical Background of Oblateness

Oblateness refers to the amount by which the equatorial radius of a planet
differs from its polar radius. A common cause of planet shape distortions are
tidal forces which induce a bulge towards the host star (Love, 1911). In our
case we are interested in oblateness induced by a planets rotation, which
causes its equator to bulge out perpendicular to its rotation axis. Standard
notation uses the parameter $f$ for oblateness, which is defined as:

$f=\frac{R_{eq}-R_{pol}}{R_{eq}}$ (1)

where $R_{eq}$ and $R_{pol}$ are the equatorial and polar radii of a planet
(Murray & Dermott, 1999). Thus oblateness ranges from 0 (for a perfectly
spherical planet) to 1 (for a planet completely flattened along the radial
direction). In principle $f$ could also be negative for the case of a prolate
planet, whose polar radius is larger than its equatorial radius. While we
restrict the oblateness to be a positive value in this work, we do allow its
rotational axis to be tilted up to 90 degrees relative to its orbital axis,
which would cause an oblate planet to appear prolate. For the planets in our
own solar system, $f$ ranges from near zero for the rocky planets such as
mercury and earth, to approximately 0.065 and 0.1 for Jupiter ($P_{Rot}=9.93\
\mathrm{hrs}$) and Saturn ($P_{Rot}=10.65\ \mathrm{hrs}$)
respectively111Values taken from
https://nssdc.gsfc.nasa.gov/planetary/planetfact.html. The rotational period
of a planet and its oblateness are related by (Hubbard, 1984)

$P_{rot}=2\pi\sqrt{\frac{R_{eq}^{3}}{GM_{p}\left(2f-3J_{2}\right)}}$ (2)

where $M_{p}$ is the planet’s mass, and $J_{2}$ is its quadrupole moment,
which is around 0.015 for Jupiter and Saturn, around 0.001 for Earth and Mars,
and $<$ 0.00001 for Mercury and Venus.

The planet may also be tilted relative to its orbital plane by two obliquity
angles $\theta$ and $\phi$. When an ellipsoid is rotated in three dimensions
and projected onto the sky plane we observe an ellipse which we can
characterise by a projected oblateness $f_{\mathrm{proj}}$and a projected
obliquity $\theta_{\mathrm{proj}}$, as shown in the left panel of figure 1.
The projected values are related to the true values by

$f_{\mathrm{proj}}=1-\sqrt{\mathrm{sin}^{2}\theta^{\prime}+(1-f)^{2}\mathrm{cos}^{2}\theta^{\prime}}$
(3a) $\mathrm{tan}\theta_{\mathrm{proj}}=\mathrm{tan}\theta\mathrm{sin}\phi$
(3b)
$\mathrm{cos}^{2}\theta^{\prime}=\mathrm{sin}^{2}\theta\mathrm{sin}^{2}\phi+\mathrm{cos}^{2}\theta$
(3c)

(a) The number of confirmed transiting planets whose lightcurve would vary by
a given amount if they had a certain level of oblateness, with a one minute
integration time. The x-axis represents the Root-Mean-Square (RMS) deviations
during ingress. The solid, dashed, dotted, and dash-dot lines represent
obliquities of 0, 45, 60, and 90 degrees respectively.

(b) For the same sample of planets as in a), we keep the oblateness fixed to
0.2 in all cases. The curves correspond to different observing cadences which
smooth out the oblateness signal.

Figure 2: The expected oblateness of confirmed transiting planets

This implies that for most planets what will be observed is an oblateness
value that is smaller than their true oblateness, diminishing the deviations
from a spherical planet. If one assumes that the orbital axis is randomly
oriented (i.e. uniformly distributed on the unit sphere), then the median
projected oblateness would be between 75% to 70% of the true oblateness (see
appendix A). Constraints on the amount by which a planet’s rotation can be
tilted away from the normal vector of its orbital plane would allow us to put
limits on the effect of projection. However in our own solar system we observe
Uranus to have its rotational axis oriented by $97.77^{\circ}$ relative to its
orbital plane, providing us with the two extremes of perfect alignment (as is
the case for Mercury) as well as perfect misalignment, along with several
intermediate values.

The observable effect of an oblate planet is illustrated in the right panel of
figure 1. During ingress and egress of a transit, two planets with differing
shapes will block differing amounts of light from their host star. More
specifically, the time-dependence of the amount of light they block will vary,
which affects the shape of the lightcurve as the planet enters and exists the
limb of the stellar disk. Once the planet is fully inside the disk of the
star, the dominating factor in the amount of light blocked is the total
projected area of the planet. Stellar heterogeneities such as limb darkening
and star spots would still produce differing signals for spherical and non-
spherical planets, however the effect is much weaker when compared to the
differences during ingress and egress, which may be be upwards of several
100’s of parts per million (ppm). It is thus critical when searching for
oblateness to ensure that the limbs of the transit are sufficiently sampled in
time.

###  2.1 The observability of oblateness for known planets

Currently, the detection of oblateness is limited to a handful of tentative
measurements among the most suitable planet candidates. Previous studies such
as Zhu et al. (2014) and Biersteker & Schlichting (2017) have used the Kepler
survey as their pool of observations, which have provided us with some of the
highest quality, as well as longest baseline, transit observations to date.

The goal of this section is to quantify the average strength of an oblateness
signal that one expects to see across all known planets. This allows us to
determine both the number of expected planets whose oblateness may be measured
with upcoming instruments, as well as the potential to place upper limits on
the oblateness of planets if no such signal is seen at increased signal to
noise ratios (SNR). The amplitude of an oblateness signal depends on
properties such as the planet-to-star area ratio, inclination, impact
parameter etc. in non-linear ways, and so planets must be assessed on an
individual basis. We compiled a list of all known transiting planets, taken
from the exoplanet archive table of confirmed planets (NASA Exoplanet Archive,
2019). As will be discussed in section 5.3, it is expected that planets with
short orbital periods ($\lesssim$ 15 days) will be rotating too slowly to be
oblate. We thus cut out planets with a period of less than 10 days, in order
to exclude planets which we confidently do not expect to exhibit oblateness so
as to not inflate the number of potential candidates. We also exclude planets
which do not have reported values for any of the transit parameters necessary
to calculate a lightcurve model.

This leaves us with a sample of $\sim$900 planets. For each of these we then
generate synthetic lightcurves for a spherical planet (using the Batman
package of Kreidberg (2015)), as well for varying levels of oblateness. The
oblateness is kept to be less than 0.5 ($r_{eq}\leq 2r_{pol}$) for stability
considerations. Equating the centripetal and gravitational forces at the
planets surface

$\frac{GM}{r_{eq}^{2}}=\frac{v^{2}}{r}$ (4)

and substituting equation 2 into $v=2\pi r_{eq}/P_{rot}$ to relate the
velocity at the surface to the oblateness parameter f, we find

$f_{crit}=(1+3J_{2})/2$ (5)

which is approximately 0.5, for $3J_{2}<<1$. To obtain values of f larger than
0.5, the necessary rotational velocity of the planet would overcome its
gravitational attraction and become unstable.

We generate 1-minute cadence lightcurves of oblate planets using a Monte-Carlo
integration routine, based off appendix A of Carter & Winn (2010a), which
includes a quadratic limb darkening law222Using the corrected versions of
equations (B5) and (B6) noted in Zhu et al. (2014),
$r=\sqrt{u+(1-u)(a_{1}/a_{2})^{2}}$ and $\theta=(1-v)\theta_{1}+v\theta_{2}$.
The uncertainty in the model itself is at the level of 1-2 ppm (set by the
number of points used in the integration).

When generating oblate lightcurves of a given planet, we fix the transit
parameters to their measured values. The only parameter we vary is the
equatorial radius of the planet, which we adjust so that the depth of the
transit remains constant. The transit depth of a planet of oblateness f is
given by

$\delta_{ob}=r_{eq}r_{polar}=r_{eq}r_{eq}(1-f)$ (6)

and so if we were to keep the radius the same, the depth between an oblate and
spherical planet would differ in proportion to f. Given that during the center
of transit there is little to no information regarding the shape of the planet
(depending on the degree of limb darkening) we set
$\delta_{ob}=\delta_{circ}$, which gives the relation between the equatorial
radius of an oblate and spherical planet as

$r_{ob}=r_{circ}/\sqrt{1-f}$ (7)

We calculate the difference between the spherical and oblate lightcurves, and
then calculate the root mean square of the difference during ingress and
egress, which we show in figure 2(a). Alternatively, we could have chosen the
amplitude or maximum difference between the lightcurves. However for many
planets which have low time-sampling during ingress / egress, it would be
misleading to report what is usually a single high point as the overall
strength of the oblateness signal.

We find that between 10ppm - 100ppm is when most planets begin to exhibit
oblateness. For nearly maximal value of f$\sim$0.5 we find as many as 50
planets with RMS deviations above 100ppm. However this is an extreme case, and
as mentioned due to projection effects it is unlikely many planets would
display such a high level of oblateness. At a signal amplitude of
approximately 50ppm, there are on the order of 30 planets which would exhibit
variations at this level if they were to have an oblateness factor of 0.2, and
almost double this number of planets if they had an oblateness of 0.3. Going
down to 10ppm, we see these numbers roughly double again.

In addition to oblateness, we also vary the obliquity angle of the planets. We
find that when considering a bulk population the obliquity does not
significantly change the number of observable planets at a given signal level.
As shown in the right panel of figure 1, it is possible for the signal to vary
significantly between ingress and egress, indicating that transit asymmetry
may be a potential sign of oblateness. If we calculate the RMS separately for
ingress and egress and compare them, we find that the ratio of their values is
on average between 0.8 and 0.9, depending on the oblateness value. For planets
with sufficiently sampled limbs this effect would be an additional marker of
oblateness.

We also study the effect of observing cadence (i.e. the integration time
between data points), which is shown in figure 2(b). Given that oblateness
primarily causes variations during ingress and egress, the number of
observations during these parts of the transit directly impacts the ability to
detect oblateness. We hold the oblateness level fixed at 0.2 and alter the
duration of the observing cadence, ranging from 1 minute to 30 minutes, as a
reflection of the Kepler observing cadence modes. We see in this case that for
a given signal level, the number of planets which exhibit oblateness decreases
by more than half as the integration time increases. We note that for most of
the planets in our sample, a one minute observing cadence corresponds to about
5 - 50 data points during the transit limbs, with a median of 16 data points.
The effect of time sampling, which smooths and distorts the shape of a
transit, is often not considered when discussing oblateness. The significance
of this is analysed further in section 4.2.

###  2.2 Where are all the oblate planets?

We have shown that as many as several tens of exoplanets should have levels of
oblateness in the 30-100ppm range at observing cadences of $<$ 10 minutes.
While a small signal, this is certainly within the capabilities of facilities
such as Kepler, as shown in figure 5 of Koch et al. (2010). To date, there is
only a small handful of planets which have had tentative detections of
oblateness, as mentioned in section 1. This begs the question of why there
haven’t been a larger number of conclusive oblateness measurements. One
explanation for this is that exoplanets are simply not oblate, which would
make the case of Saturn’s oblateness an extreme outlier, or that oblate
planets happen to be found around dimmer stars with lower photometric
precision. Another option which we will explore in the next section is that
the transit parameters of associated with a spherical planet model are able to
mimic the signature of oblateness and thus compensate for it.

##  3 Mapping the Degeneracies Between Oblateness and Transit Parameters

As shown in the previous section, the typical deviations induced by oblateness
will be in the best case a few hundred ppm, and much more typically well below
100 ppm. The additional fact of the signal being confined solely to the limbs
of a transit means that for many planets, depending on the duration of the
transit and the cadence at which the observations are taken, only a handful of
data points may be useful in detecting oblateness. It is thus not surprising
that when analysing an exoplanet transit the standard practice is to assume a
perfectly spherical planet. Nevertheless, we propose that even in cases of
non-detection oblateness may have a significant effect on parameter retrieval.

Figure 3: Corner plot showing the deviations for each parameter from their
truth values, in units of standard deviations of the posterior. Colors
represent different levels of oblateness. The black circles show boundaries of
1, 2, and 3 standard deviations.  Figure 4: Results of fitting a spherical
planet model to simulated data of an oblate planet, using orbital parameters
for the 100 largest (in terms of planet to stellar radius) Kepler planets with
periods greater than 10 days. The left vertical axis in each plot represents
the number of standard deviations between the truth and fitted value for the
planet to stellar radius, impact parameter, inclination, quadratic limb
darkening parameters, eccentricity parameters, and transit center. Dark and
light regions represent where 65% and 95% of the planets in the sample appear.
The right vertical axis and solid red line show the percentage of planets for
which the retrieved parameters deviate by $>$2 standard deviations. The bottom
right panel shows the median values of the reduced $\chi^{2}$ in red, as well
as the median Bhattacharyya distance in blue.

###  3.1 Bias on retrieved parameters during a transit analysis

Parameters of a model are often degenerate with one another, and given the low
signal strength of oblateness it’s possible that other transit parameters may
be able to compensate and mask an oblateness signal at current SNR. Assuming a
planet to be perfectly spherical could then introduce a bias when fitting for
parameters such as inclination or the semi-major axis of the planets orbit.
The goal of this section is to quantify such an effect through a series of
simulation-retrieval studies. We note that in Barnes & Fortney (2003) a
similar analysis was done to study the biases of assuming an oblate planet to
be spherical. However in that case the effect was only studied for a single
planet, HD209458b (Charbonneau et al., 2000), whereas our aim is to expand
this analysis to a wider set of planets to understand the effect more
globally. We follow the methodology of de Wit et al. (2012) which looked at
the degeneracies between the orbital properties of a planet and its potential
to have non-uniform surface brightness, measured through minor variations
during eclipse.

We first generate a set of lightcurves for a sample of oblate planets. The
planet parameters we use come from the exoplanet archive list of confirmed
transiting planets, as was done in the previous section. The parameter space
of possible planet and orbital configurations is large, and so we restrict our
study to known planets, which at the time of writing is a sample of greater
than 4000 objects, and from these we take the 100 planets with the largest
planet to stellar radius ratios (since the oblateness signal scales
proportionally with radius). For those planets, we generate 1 minute cadence
lightcurves with oblateness value from the range of 0, 0.1, 0.2, 0.3, and 0.5.
We add Gaussian noise with a standard deviation of 50ppm for each data point.
We then run a Monte-Carlo Markov-Chain (MCMC) fit to the synthetic data using
the emcee package (Foreman-Mackey et al., 2012), assuming a spherical planet
model during the fit.

###  3.2 Description of Fit

We choose to fit for planet-to-stellar radius ratio, transit impact parameter,
inclination angle, eccentricity, argument of periastron, as well as two
parameters of a quadratic limb darkening model, using the q1 & q2
parametrization of Kipping (2013) and finally center of transit. We hold the
orbital period fixed. We note that some of these parameters, specifically the
eccentricity, could be constrained by radial velocity observations which would
be unaffected by planetary shape variations. Additionally, for the fits
described in this section we hold the obliquity angle of the planet fixed at
45 degrees. We ran a similar analysis for a fixed obliquity of zero which we
discuss further down.

Regarding the radius of the planet, we do not fit for it directly but rather
scale it by the appropriate oblateness factor as previously mentioned. If we
fit for the radius directly, we would see deviations on the order of 20% or
more as oblateness increased (in order to match the transit depth), however
this would be somewhat misleading. While it is true that the equatorial radius
would differ by such an amount, the surface area (and thus total occulted
stellar area) would remain relatively unchanged. This means that a 20% change
in equatorial radius between a spherical and oblate planet does not represent
the same magnitude of change between two spherical planets whose radii differ
by 20%. By fitting for the scaled radius, any deviation we see is a more
meaningful difference in quantities such as the total surface area or volume
of the planet.

For a given planet we compare the posterior distributions of each fitted
parameter to the known value that was used to generate the simulated data. We
calculate the median value of the posterior distributions as well as the
standard deviation. We then calculate by how many standard deviations the
fitted parameter value is from the true parameter used when generating the
data. If one assumes gaussian uncertainties for the fitted parameters, then
for an appropriate model it would be expected that $\sim$68% of the time the
true parameter should be within one standard deviation, $\sim$95% of the time
within two standard deviations and $\sim$98% of the time within three standard
deviations.

###  3.3 Results of injection-retrieval tests

The compiled outputs of our set of fits is shown in figures 3 & 4. As
expected, for low oblateness values the fitted parameters cluster around the
true value, with a spread of $\sim$ 2 standard deviations. This is the
statistically expected result for a model which matches the data.

####  3.3.1 Biases on radius, impact parameter, and inclination

We observe significant biases with clear trends for the equatorial radius
$r_{p}/{r_{\star}}$, impact parameter (b), and inclination (incl). For
oblateness values greater than 0.2, we see that the first two of these
parameters become biased towards smaller values, while inclination is biased
towards larger values (i.e. the orbits tend towards the observer line of
sight). We also see the distribution of fit parameters widen, which could be
an indication that the errors on these parameters are being underestimated.

For ($r_{p}/{r_{\star}}$) we find that 20% of the fits deviate by $>$
2$\sigma$ at f = 0.2, with this number increasing to $>$60% as the oblateness
increases. For the impact parameter we find that 10% of the fits deviate by
$>$ 2$\sigma$ increasing to 40% at f = 0.5, and for inclination we find that
from f = 0.2 to f = 0.5 the percentage of fits $>$2$\sigma$ rises from
$\sim$5% to 40%.

A potential explanation for this is the effect of oblateness at the very
beginning and end of transits. Compared to a spherical planet, and based on
its orientation, an oblate planet will begin transiting slightly earlier and
end transiting slightly later. This leads to a change in the overall duration
of transit, which we see here is being compensated for by the three parameters
which directly affect the length of the transit chord.

In figure 3 we note a significant correlation in particular between the impact
parameter and inclination of the planet. We observe a similarly strong
correlation between the inclination and radius ratio.

For the specific case of HD209458b, it was found in Barnes et al. (2007) that
when fitting a spherical planet model to a simulated lightcurve of an oblate
planet the impact parameter tended towards a critical value of $b=0.707$ in
order to account for the change in ingress/egress duration. We do not observe
a similar trend for the larger sample of planets and parameters we have
studied in this work. One important distinction between the two analyses is
that we have included additional parameters in our fit, namely the
inclination. The effect of this is that we have additional ways to account for
the change in transit durations. Thus we find differing trends in the
combination of impact parameter and inclination which can account for the
variations between spherical and oblate lightcurve models.

####  3.3.2 Biases on limb darkening coefficients

For both limb darkening parameters q1 and q2 there is no clear upwards or
downwards trend, although we do see a widening of distributions (although with
a smaller amplitude of deviations compared to the other parameters). If a
change in transit duration is indeed the cause of any observed biases, this
would make sense given that the limb darkening parameters are unable to alter
the length of the transit chord and are less strongly biased. Despite a lack
of a clear trend however, we still find that a significant number of fits
deviate by more than 2 standard deviations, indicating that the limb darkening
coefficients are still affected by the change in transit shape to some degree.

For q1 we find that 10% of the fits deviate by $>$ 2$\sigma$ at f = 0.2,
increasing to 60% as the oblateness increases towards 0.5. For q2 we see a
larger initial value, with $\sim$20% of fits deviating by $>$ 2$\sigma$ at f =
0.2, increasing to $\sim$40% at f = 0.5.

####  3.3.3 Biases on transit center

In the case of non-zero obliquity, we observe a significant shift in transit
center for large oblateness values. As will be discussed in section 5.2, this
is an expected effect due to the asymmetry in ingress and egress for a tilted
planet. As previously mentioned, we ran a similar analysis for a fixed
obliquity of zero. In this case, the transit will always be perfectly
symmetric about its center, and thus no shift in transit center is expected
which is what was observed. We note here that the direction of the shift is
always the same. This is due to the fact that the obliquity is the same for
all planets, and we assign the impact parameter to positive values. Thus the
asymmetry induced in the lightcurve will always be in the same direction with
regards to ingress and egress. For f = 0.2 $\sim$ 10% of fits deviate by
$>$2$\sigma$, increasing to 40% as f approaches 0.5.

####  3.3.4 Biases on eccentricity

For eccentricity and the argument of periastron $\omega$ (parameterized as
$\sqrt{e}cos{\omega}$ and $\sqrt{e}sin{\omega}$) we find that the retrieved
values are statistically consistent with the values used to generate the data
(i.e. all samples are $<$3$\sigma$). We note however that this is likely to be
a limitation of the data sets being analysed, which are single transit
observations. Eccentricity can be more tightly constrained when additional
information such as stellar density (Dawson & Johnson, 2012) and planetary
eclipses are used (Winn et al., 2007). We thus would expect to obtain
significant biases on these parameters from perturbed ingress/egress for the
increased precision derived from joint fits of primary eclipses, secondary
eclipses, and/or radial velocity measurements–similarly to the point made in
Fig. 4 of de Wit et al. (2012).

####  3.3.5 Goodness of Fit Metrics

In addition to having measured a clear bias for certain parameter when
assuming a spherical planet model, we also check the goodness of fit statistic
to asses if, given current precisions, there exists a range of oblateness
leading to perturbations in the transit shape which cannot be compensated for
by the orbital parameters. We show in the bottom right panel of figure 4 the
median value of the reduced $\chi^{2}$ across all planets in the sample for
differing oblateness values. We note that this is measured only using data
points during ingress and egress of transit, as the rest of the lightcurve is
expected to be the same for both spherical and oblate planets. As shown, the
value remains close to one, increasing only for values of oblateness
approaching 0.5 (although still remaining below $\sim$1.5). In the same panel
we also show the mean value of the Bhattacharyya distance ($D_{B}$), which is
a metric for comparing the similarity of two distributions defined as

$D_{B}(p,q)=-\mathrm{ln}\left(\Sigma_{i}\sqrt{p(\theta_{i})q(\theta_{i})}\right)$
(8)

for two probability distribution p and q over a set of parameters $\theta$.
The value shown for a given oblateness is a comparison between the multi-
dimensional posterior distribution of that value and the posterior
distribution for zero oblateness. We see a clear trend to higher values as the
oblateness increases, indicating a statistically significant variation between
the two distributions. Thus we find that assuming planets to be spherical can
statistically bias results in a way that does not set off the most common
alarm for improper models.

##  4 The sensitivity of oblateness retrieval & planet populations

Figure 5: The results of MCMC retrievals to two populations of planets, using
a model which accounts for oblateness and obliquity. On the left is a
population where each planet has an oblateness value of 0.2, and the the right
plot is a population of planets with zero oblateness. For each group, we test
difference levels of observing cadence (dt in minutes) and Gaussian noise
level ($\sigma$ in ppm)

Measuring the oblateness of an individual planet is difficult and requires
very precise and well-sampled observations. The approach taken in this section
is to study ensembles of planets in order to determine if planets can be
separated into populations of oblate and non-oblate planets. For example, as
given by equation 12, the tidal spin-down timescale varies strongly with
orbital period (to the $\mathrm{4^{th}}$ power). Planets below a period of 15
days are expected to have very short spin-down timescales, implying they will
be tidally locked to their host star and will have slowed down enough to have
effectively zero oblateness. For planets with an orbital period above 15 days
with much longer spin-down times, they are expected to have retained a high
rate of rotation from their formation and thus have non-zero oblateness. In
addition to searching for oblateness across a population, we also measure the
sensitivity of oblateness fits to the quality of the data, which we
characterise by SNR and time sampling of data points.

###  4.1 Synthetic populations

We generate 200 planets whose periods, radii, semi-major axis, and
inclinations are drawn from the population of known transiting planets. For a
given choice of observing cadence, oblateness factor, and noise level we
simulate an oblate lightcurve as in previous sections. We do this process for
two populations of planets with identical planet parameters, with one having
zero oblateness and the other having an oblateness of 0.2 for all planets.

Once we have these lightcurves, we run an MCMC retrieval on them using the
Allesfitter package (Günther & Daylan, 2021, 2019). We include in Allesfitter
our own oblateness model, which allows us to fit for the oblateness parameters
f (oblateness) and $\theta$ (the obliquity of the planet).

We show in figure 5 the results of this analysis, in particular the posterior
distributions retrieved for the oblateness. We note that in our analysis we
reparametrize f and $\theta$ as $\sqrt{f}\mathrm{cos}\theta$ and
$\sqrt{f}\mathrm{sin}\theta$, similar to the parametrization often used when
fitting for eccentricities (Van Eylen & Albrecht, 2015). These reparametrized
versions of the shape parameters do not have hard boundaries at 0 and instead
vary from -1 to 1, which removes the bias of forcing oblateness to be
positive.

We first note that in the case of zero oblateness (the right panel), the
general trend is a flat posterior below a certain value, in this case around
$f=0.25$, followed by a sharp decrease beyond that. The model is able to
confidently rule out high oblateness values, which in turn can be used to rule
out large rotational periods of planets. Below the cutoff point the model is
less sensitive to oblateness. We see that for the most constraining set of
parameters, namely a short observing cadence of dt = 1 minute and a small
noise level of 50 ppm, the model is able to exclude smaller oblateness values
more confidently.

In the left panel of figure 5 we see two cases emerge, depending on the noise
properties and time sampling. For the shortest time sampling and lowest noise
level, the model is able to retrieve the population oblateness level of 0.2.
As the noise level increase along the with the time sampling, the posterior
gradually shifts to that of a zero oblateness population. Thus we see that
between 50-100 ppm, and between 1min - 10 min observing cadence there is a
cutoff point where oblateness shifts from being confidently detectable to
completely undetectable.

###  4.2 Sensitivity of oblateness fits to data quality

In order to investigate the shift in retrieval seen, we sample a grid of data
quality, which we characterise by two metrics. One metric is the signal-to-
noise ratio, which in our case we take to be the transit depth divided by the
amplitude of Gaussian noise. The other value we use to quantify the data
quality is the number of data points observed during ingress, which is
calculated as the duration of transit divided by the observing cadence. For
each data point sampled, we run an MCMC retrieval as described above, and
calculate the oblateness value which encompasses 95$\%$ of the posterior
distribution (which is often asymmetric due to the boundary at f = 0).

Figure 6: The median 2$\sigma$ oblateness upper limit values across a sample
of planets retrieved when fitting an oblate planet model to synthetic data
sets with zero oblateness (values shown in grid cells). The x axis indicates
the transit depth divided by the amplitude of Gaussian noise. The y axis shows
the duration of transit ingress divided by the observing cadence.

The results of this are shown in figure 6. In this case the data being fit has
an oblateness value of zero, i.e. spherical planets. Thus we are demonstrating
the ability of the model to confidently rule out non-zero oblateness. We see
that for low SNR values and a low number of transit data points, the model is
unable to rule out values of oblateness below 0.3 at the one standard
deviation level, and f = 0.4 at the 2$\sigma$ level. It is only for very high
SNR values $>$ several hundred and more than at least 10 transit points that
the model is even able to rule out oblateness values of 0.2 or greater. We
note that the oblateness value within each cell varies to a certain degree
amongst the different planets. This is due to the fact that properties such as
the impact parameter also play a role in measuring oblateness, however here we
focus on the bulk results marginalized across such parameters.

The significance of this is that the majority of transiting planet
observations are completely degenerate between a spherical planet and a planet
having an oblateness of $\lesssim$ 0.25. In the context of section 5.1 for
example, this implies that analysis which rely on relative density errors
below $\sim$ 15$\%$ are over-confident, in the case of planets which have the
potential to be oblate (i.e. planets which are unlikely to be tidally locked).

###  4.3 Oblateness fit of long and short period Kepler planets

Figure 7: Oblateness posteriors when fitting confirmed Kepler planets,
separated by orbital periods.

We now turn to the potential oblateness of planets discovered by the Kepler
survey. The reason for picking this as our sample is the large number of
planets it contains, the fact that many of them were observed at a short (1
minute) cadence, the duration of the survey which provides a large number of
transits for most planets, and finally the SNR it provides for most planets,
which far surpasses ground based and many spaced-based transit observations.

In contrast to previous studies which have searched for oblateness in
individual planets, our goal is to see if it is possible to highlight
populations of planets which globally either do or do not exhibit oblateness.
Given that planetary rotation is the root of the signal, we propose to
separate planets into fast and slow rotators, which would then translate into
populations of high and low oblateness planets.

We download the short cadence lightcurves of all planets with radii larger
than a cutoff value of 6 $R_{e}$, to exclude rocky planets with low tidal
dissipation factors (and thus extremely short circularization timescales). We
do this using the lightkurve package (Lightkurve Collaboration et al., 2018),
and extract the PDCSAP lightcurves. We then run them through a suite of post-
processing and normalizing steps. We first remove points which have bad data
quality labels. We next identify all transits of a system, including those
with multiple planets. We chose not to exclude these systems, although we do
remove observations where two transits come within a multiple of 3 transit
durations of one another. We then apply a Savitzky-Golay filter to remove long
term variations in the lightcurve. We mask the known transits while doing this
in order to not have them bias the detrending. We next correct for systems
which exhibit large TTVs, which in some cases by on the order of tens of
minutes. We do by this running an MCMC fit for each individual transit,
holding fixed all known transit parameters and allowing only the transit
center and depth to vary.

Once we have measured transit centers for all planets, we repeat the previous
detrending steps, using the shifted transit centers to accurately mask out the
transits when fitting for low frequency lightcurve trends. Once this is done
we phase fold the data into a single lightcurve, which is what is finally used
in the fitting routine. We run an MCMC fit to the data utilizing the
Allesfitter package as we did in the previous sections.

We show in figure 7 the stacked posterior distributions for planets with
periods either above or below 15 days. We find similar distributions as in the
previous section for the cases in which the fit was unable to identify
significant oblateness. The fact that both distributions look similar indicate
that the sample of Kepler planets does not exhibit an oblateness bias or
grouping based on their period. We analyzed different groups of periods,
cutting below 10 days or above 20, and similarly found no biases or
differences between any groups.

We also see the same cutoff at oblateness values above $\sim$0.25, again as in
the previous section. This indicates that the MCMC disfavors such high levels
of oblateness. High oblateness values in general are unlikely, due to the
projection effects described in section 2 and Appendix A. In section 2.1 we
showed the oblateness signal induced among known planets, and show that very
few planets produce a signal beyond 100ppm. However for many planets in the
Kepler survey the noise level is on the order of 100s of ppm. Thus while we
are able to exclude high values of oblateness in general, we find that the
sensitivity factors described in the previous subsection limit the ability to
detect or rule out more moderate values of oblateness.

We do note a slight increase in the posterior for oblateness values $<$ 0.05
for planets with orbital period $>$ 15 days. This is again most likely due to
increased sensitivity discussed in the previous section, given that planets
with larger orbital periods will have increased sampling during ingress /
egress and also tend to be larger. Thus we do not attribute this as indicating
that larger period planets have lower oblateness.

##  5 discussion

So far we have studied the effects on standard transit analysis when allowing
for an exoplanet to be oblate instead of spherical. If not properly taken into
consideration, ignoring the potential of a planet to be oblate could introduce
additional biases into more sophisticated exoplanet analysis. This highlights
a main point of this work, which is that as the quality of data continues to
improve, features such as oblateness which were once undetectable may begin to
bias and effect a wider range of conclusions one may make about explanatory
systems. In this section we consider the measurement of bulk densities of
planets, as well as transit timing variations and the ways in which oblateness
may impact any conclusions drawn. We also describe a metric for quantifying
which planets are more or less likely to exhibit oblateness, based not only
signal strength but also on their likelihood to be oblate in the first place.

###  5.1 The relative uncertainty of bulk densities

When discussing radii in section 3, care was taken when comparing a spherical
planet described by just one radius to an elliptical planet with two radii.
However for a quantity such as the density of the planet, there is less of a
distinction since the bulk density is typically the quantity of interest. We
will illustrate that any analysis which utilizes precise density measurements
may be subject to biases due to the shape of the planet.

Consider a planet which has an oblateness f, which induces oblateness at a
level that is too low to be detectable or even to induce the biases previously
mentioned. The depth - radius relation for an oblate planet is given by
$\delta=r_{eq}r_{pol}=r_{eq}^{2}(1-f)$. This would translate to a measured
planet volume of

$Vol_{ob}=\frac{4\pi}{3}r_{eq}^{2}r_{pol}=\frac{4\pi}{3}r_{eq}^{3}(1-f)=\frac{4\pi}{3}\frac{\delta^{3/2}}{(1-f)^{1/2}}$
(9)

Assuming mass measured by radial velocities to be independent of shape, we can
compare the inferred volumes of a spherical and oblate planet, and thus their
densities, as

$\frac{\rho_{circ}}{\rho_{ob}}=\frac{4\pi}{3}\frac{\delta^{3/2}}{(1-f)^{1/2}}\frac{3}{4\pi}\frac{1}{\delta^{3/2}}=\frac{1}{(1-f)^{1/2}}$
(10)

For an oblateness of even just $f=0.1$, this would translate to a true density
that is 5% smaller than the density inferred for a spherical planet. An
oblateness of $f=0.2$ would mean a discrepancy of almost 11%.

Figure 8: The cumulative number of planets with a fractional density error
below a certain value (blue line). The black lines show the percentage by
which the density would change if a previously assumed spherical planet were
to be oblate to a certain degree, calculated using equation 10. The red line
indicating a value of f = 0.25 is in reference the limits of current
observations discussed in section 4.2.

In figure 8 we compare this oblateness induced density variation to the
relative uncertainty of the density of known planets. We take all planets from
the exoplanet archive which have a reported value of both physical radius and
mass, as well as reported uncertainty values. We calculate their relative
density errors, and find for an oblateness of 0.2, up to 20 planets would have
their density differ by more than their relative uncertainty. For an
oblateness of 0.3, that number rises to almost 60. It was demonstrated in
section 4.2 that at low data sampling and/or transit SNR it is impossible to
rule out oblateness values less than 0.3 at the 2$\sigma$ level. This in turn
imposes an accuracy wall on relative uncertainty of a planet’s bulk density at
the level of 10$\sim$15%.

###  5.2 TTV-like Signal induced by Oblateness

(a) The distribution of retrieved TTV amplitudes induced by asymmetry stemming
from oblate planets with f = 0.2 and $\theta=\pi/4$. The y axis shows the
impact parameter of the transiting planets, which we expect to be correlated
with the asymmetry signal.

(b) The orbital periods of a companion that would be needed to induce the TTV
signals seen, based off the precession period and possible orbital resonance.

Figure 9: The oblateness induced TTV effect

Thus far we have considered the oblateness of a planet to be constant, however
it is possible that a planet with a non-zero obliquity angle will precess over
time as it orbits its host star and exchanges angular momentum. For uniform
precession in a fixed orbit, the precession period is given by (Ward, 1975)

$P_{prec}=13.3\ \mathrm{years}\times\left(\frac{P_{orb}}{15\
\mathrm{days}}\right)^{2}\left(\frac{10\mathrm{hr}}{P_{\mathrm{rot}}}\right)\left(\frac{\lambda/J_{2}}{13.5}\right)\frac{1}{\mathrm{cos}\theta}$
(11)

where the value of 13.5 comes from the estimate of $\lambda/J_{2}$ for Saturn
(Ward & Hamilton, 2004). In this equation $\theta$ is the obliquity angle of
the planet, which we recall is the angle between the polar axis of the planet
and its orbital plane. As mentioned in section §2, as the orientation of an
oblate planet changes in space, an observer will see different elliptical
projections on the sky plane. This will lead to a time varying projected
oblateness and obliquity, which will in turn lead to a time varying transit
shape.

The focus of Biersteker & Schlichting (2017) was to try and detect oblateness
due to do variations in the depth of a planets transit, caused by a time
varying projected surface area (including a change in both the projected
oblateness as well as equatorial radius). A secondary effect which was first
mentioned in Carter & Winn (2010b) is that the duration of the transit may
change as well, which could mimic the appearance of a transit timing variation
(TTV) signal.

If an oblate planet has a non-zero impact parameter and obliquity, then it
will exhibit asymmetry in ingress and egress (as in figure 1). Thus when
attempting to fit a spherical transit model to the lightcurve of an oblate
planet, the result could be a perceived shift in the center of the transit.
This effect is typically attributed to another planet in the same system
gravitationally perturbing the observed planets orbit (Agol & Fabrycky, 2018).
Thus attempting to fit a spherical planet model to a precessing oblate planet
could potentially lead one to conclude the existence of an additional planet
causing timing variations.

We perform a set of injection retrieval tests, attempting to measure a shift
in the center of transit. We simulate a transit lightcurve of an oblate planet
across multiple epochs while altering the obliquity over time according to its
precession period. As the true obliquity of the planet changes over time, we
also update the projected obliquity and projected oblateness. We must
additionally account for the change in projected equatorial radius of the
planet, which we do using equation 4 of Carter & Winn (2010b). The formula
provides the relative transit depth as a function of oblateness and
inclination parameters.

We use the same sample of planets as in previous sections in order to
determine the degree to which the current population of planets may be
affected. We note that the precession period given in equation 11 is only an
estimate, and further effects such as interactions with other planets in the
same system could disrupt or alter the precession in someway (Saillenfest et
al., 2019). We hold the parameter ${\lambda/J_{2}}$ fixed to the value of 13.5
for Saturn, and note that the precession period scales linearly with it. We
sample the obliquity of the planet from three value spaced out between 0 and
PI/2. In the case of of either zero or maximal (PI/2) obliquity, there would
be no ingress/egress asymmetry and thus no induced TTV signal. For these
values we find that most precession periods range from 10 - 100 years,
although there are some planets with potential periods of only a couple years
or less. For most planets such an effect will thus occur on too long of a
timescale to be observable. As the field of exoplanet enters its $4^{th}$
decade however, it becomes more and more feasible to observe such a long
baseline transit shift and falsely attribute it to a long period companion.

Many of the trials return a constant transit center as the best fit,
indicating no detectable variation. For a significant number of cases however
the best fit is a clear sinusoidal variation in the transit center, as one
would expect from TTVs. For these cases we fit a sinusoid to the time varying
signal and retrieve the amplitude of the variation. In figure 9(a) we show the
distribution of these amplitudes, in units of seconds, as a function of the
impact parameter of the transit. We note a correlation with the impact
parameter, which is extremely tight for values less than 0.15, and then widens
above that. We find that for an underlying oblateness of 0.2 and an obliquity
of 45 degrees, the scale of the deviations is on the order of several seconds,
with the peak value being upwards of 15 seconds (and thus 30 seconds peak to
peak). While small, transit deviations have been measured at the level of
seconds, for example in Wasp-12 b (Maciejewski et al., 2016).

In figure 9(b) we show the possible orbital periods of the companion planets
that would be required to induce the observed TTV signal. This is calculated
using equations 6 and 7 of Lithwick et al. (2012), where we substitute the
precession period of the oblate planet for the super-period of the TTV signal.
We calculate a companion period under the assumption of a 2:1, 3:2 and 4:3
period resonance. The periods we find span the range of 10 - 1000 days.

As mentioned we repeated the above analysis while varying the obliquity angle
of the planet. This was done to ensure the results are consistent with the
expectation that for obliquity angles further from 45 degrees the induced
amplitudes should be smaller due to the reduced asymmetry between ingress and
ingress. This is indeed what was found, and in particular for an obliquity of
zero/90 degrees (in which the entire transit is perfectly. symmetry about its
midpoint) there were no statistically significant transit center variations.

###  5.3 A metric to identify prime oblateness candidates

In this section we define an oblateness observability metric to determine
which planets are the most ideal for a targeted study of oblateness. It is
built in such a way to favor planets which would not only have a large
oblateness signal, but which are also likely to be oblate in the first place.
The key factor we use to determine if a planet is likely to be oblate or not
is the age of its host star. As a planets rotation slows down over time due to
tidal de-spinning, its oblateness similarly decreases. The timescale for a
planet to spin down and become tidally locked can be estimated by:

$\begin{split}\tau_{spin}=1.22\
\mathrm{Gyr}\times\left(\frac{M_{p}}{M_{J}}\right)\left(\frac{Q_{p}}{10^{6.5}}\right)\left(\frac{\lambda}{0.25}\right)\\\
\times\left(\frac{P_{\mathrm{orb}}}{15\
\mathrm{days}}\right)^{4}\left(\frac{R_{J}}{R_{e}q}\right)^{3}\left(\frac{10hr}{P_{\mathrm{rot,i}}}-\frac{10hr}{P_{\mathrm{rot}}}\right)\end{split}$
(12)

which is derived by integrating the spin down rate of a planet as described in
Goldreich & Soter (1966). The quantities $Q_{p}$ and $\lambda$ are the tidal
dissipation factor and normalized moment of inertia
$\lambda=I/(M_{p}R_{p}^{2})$. If the mass of the planet is known, one can
additionally directly calculated the rotational period required for various
levels of oblateness.

Given the uncertainty and difficulty in calculating Q, we conservatively set
it to $10^{6.5}$, which is above the estimate for the solar system values of
Saturn and Jupiter, which are typically on the order of $10^{4}$ and $10^{6}$
respectively (Storch & Lai, 2013). By over-estimating Q, we in turn will
likely be over-estimating the tidal spin-down timescale. Thus, we can be
confident that systems which are older still than this estimate will have
already spun-down and are unlikely to be oblate. For the parameter $\lambda$
we adopt the commonly used value of 0.25 for gas giants. As a comparison, the
Darwin-Radau equation333$\frac{J_{2}}{f}=-3/10+5/2\lambda-15/8\lambda^{2}$
gives a value of $\lambda=0.22$ for Saturn, and $\lambda=0.27$ is calculated
for Jupiter from Juno gravity measurements (Ni, Dongdong, 2018).

Conversely, planets whose associated star has an age that is well below 1 Gyr
can be assumed to not have had enough time to spin down and may still have a
high rotation gained during formation. We incorporate this information using a
sigmoid function, defined as $\sigma(x)=1/(1+e^{-x})$ which is a smoothed out
step function. We scale the step width by the uncertainty of the stellar age,
to allow for flexibility for spin down estimates near the stellar age.

We multiply the factor for the timescale by both the RMS variation induced by
oblateness (as described in section 2.1), and additionally by a term which is
proportional to the expected SNR of a detection, for a given waveband. This
takes the form

$\begin{split}Metric=RMS_{ob}\times
10^{-M/5}\times\sigma((Age_{\star}-\tau_{spin}^{upper})/\sigma_{Age_{\star}})\end{split}$
(13)

Where M is the magnitude of the host star in a given waveband, and the
exponential term signifies the square-root SNR estimate of a given number of
photons. The $RMS_{ob}$ term is the same ingress/egress deviation measured in
the section 2.1.

Planet |  | Period
---
(days)
| Radius
---
(Rj)
K-Mag |  | Ob. Amp
---
(ppm)
TOI-1899 b | 29.02 | 1.15 | 10.51 | 110.45
TOI-1278 b | 14.48 | 1.09 | 9.74 | 57.67
Kepler-699 b | 27.81 | 1.46 | 13.67 | 104.8
CoRoT-10 b | 13.24 | 0.96 | 11.78 | 66.27
TOI-837 b | 8.32 | 0.77 | 8.93 | 29.19
Table 1: Top ranked planets according to oblateness metric

##  6 Conclusions

In this work we have analyzed the effect of planetary oblateness on the
analysis of transit lightcurves, primarily through the use of injection-
retrieval studies of spherical / oblate planets using spherical / oblate
transit lightcurve models. We first estimated the level of variation between
the lightcurves of a spherical and oblate planet which otherwise have
identical orbital and planetary parameters. We find that for most planets with
periods greater than 10 days, this variation is at the 10-100ppm level. This
result is obtained across a sample of planets which range in size, orbital
period, impact parameter, and semi-major axis, which all affect the structure
of an oblate planets lightcurve.

Given that this is a level at which instruments such as Kepler are capable of
reaching, we then studied effects which may obscure oblateness signals,
notably compensation by orbital parameters (Barnes & Fortney, 2003; de Wit et
al., 2012). In addition to SNR, we also analyzed the effects of time sampling
on the ability to robustly measure or rule out oblateness. What is found is
that ruling out Saturn like oblateness (f = 0.1) at the 2$\sigma$ confidence
level requires a sampling of $>$ 10 data points during ingress/egress and a
transit depth which is $\geq$ 100 times the level of uncorrelated Gaussian
noise. Combining both of these effects implies that for many planets there
will be a sensitivity wall which limits the ability of a retrieval method such
as MCMC to detect oblateness. When analyzing a sample of short cadence Kepler
observations, we find this limiting factor to be in the range of f = 0.15 -
0.25. A consequence of this result is a limit on the accuracy of the bulk
density which can be measured for planets which are potentially oblate. We
find that relative density uncertainties below 10$\%$ for planets with periods
$>$ than 10 days are likely to be overestimated, given the difficulty in
ruling out oblateness levels below f = 0.2. If we extend this statement to
planets which have been distorted through other processes such as tidal
bulging, then the effect may be prevalent for planets with shorter periods as
well.

In addition to density, we also find that attempting to fit spherical transit
models to data which come from oblate planets produces fits which are both
statistically consistent with regards to the reduced $\chi^{2}$, but which
also deviate significantly from the values used to generate the data. For an
oblateness approaching 0.5, we do find that the reduced $\chi^{2}$ begins to
approach 2, indicating that the model is insufficient for the data. However in
the range of f = 0.1 - 0.3, we find that it remains close to 1, while the
Bhattacharyya distance becomes $\gtrsim$ 3 indicating deviation from the true
parameter distributions.

The planet radius, inclination, impact parameter, and transit center all
demonstrate biasing trends as oblateness increases, while the two parameters
of a quadratic limb darkening law deviate symmetrically about the true value.
The eccentricity and argument of perihelion however remain unbiased and
consistent with the truth values, even for large values of oblateness. This is
expected given that we do not consider additional information such as
secondary eclipses or radial velocity observations, which would provide a much
more accurate measurement of eccentricity (in which case we would expect
biases to appear). Additionally, time variation of oblateness due to
precession is capable of mimicking a transit timing variation signal. We note
however that this variation would be over significant timescales ($>$ 1 year)
and only at the level of 10’s of seconds.

When considering the biases and degeneracies we have found for oblateness, we
note that there is the potential to extend this to other sources of shape
variation. The models of Barnes et al. (2009) describe the effect of winds on
an exoplanets atmosphere, which would distort its outer shell into a myriad of
potential shapes which may each bias a lightcurve fit in different ways.

Furthermore, atmospheric studies of exoplanet surfaces, such as those which
leverage the eclipsing of a planet by its host star, could benefit
significantly by having the properties of the planets shape be pinned down
precisely through an analysis of its transit, where the signal is often much
stronger. Alternatively, if a planet which is potentially oblate is assumed to
be spherical, this may further bias the results of any higher order analysis
which is based off an understanding of the orbital and physical properties of
a planet derived from transit observations.

##  7 Acknowledgements

DB acknowledges support from an NSERC PGS-D scholarship as well as an FRQNT
Doctoral Research Scholarship.

## References

  * Agol & Fabrycky (2018) Agol, E., & Fabrycky, D. C. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 7, doi: 10.1007/978-3-319-55333-7_7
  * Akinsanmi et al. (2020) Akinsanmi, B., Barros, S. C. C., Santos, N. C., Oshagh, M., & Serrano, L. M. 2020, MNRAS, 497, 3484, doi: 10.1093/mnras/staa2164
  * Barnes et al. (2007) Barnes, J. R., Leigh, C. J., Jones, H. R. A., et al. 2007, MNRAS, 379, 1097, doi: 10.1111/j.1365-2966.2007.11990.x
  * Barnes et al. (2009) Barnes, J. W., Cooper, C. S., Showman, A. P., & Hubbard, W. B. 2009, ApJ, 706, 877, doi: 10.1088/0004-637X/706/1/877
  * Barnes & Fortney (2003) Barnes, J. W., & Fortney, J. J. 2003, ApJ, 588, 545, doi: 10.1086/373893
  * Biersteker & Schlichting (2017) Biersteker, J., & Schlichting, H. 2017, AJ, 154, 164, doi: 10.3847/1538-3881/aa88c2
  * Carter & Winn (2010a) Carter, J. A., & Winn, J. N. 2010a, ApJ, 709, 1219, doi: 10.1088/0004-637X/709/2/1219
  * Carter & Winn (2010b) —. 2010b, ApJ, 716, 850, doi: 10.1088/0004-637X/716/1/850
  * Charbonneau et al. (2000) Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45, doi: 10.1086/312457
  * Dawson & Johnson (2012) Dawson, R. I., & Johnson, J. A. 2012, ApJ, 756, 122, doi: 10.1088/0004-637X/756/2/122
  * de Wit et al. (2012) de Wit, J., Gillon, M., Demory, B.-O., & Seager, S. 2012, A&A, 548, A128, doi: 10.1051/0004-6361/201219060
  * Foreman-Mackey et al. (2012) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2012, ArXiv e-prints.  https://arxiv.org/abs/1202.3665
  * Goldreich & Soter (1966) Goldreich, P., & Soter, S. 1966, Icarus, 5, 375, doi: 10.1016/0019-1035(66)90051-0
  * Günther & Daylan (2019) Günther, M. N., & Daylan, T. 2019, Allesfitter: Flexible Star and Exoplanet Inference From Photometry and Radial Velocity, Astrophysics Source Code Library.  http://ascl.net/1903.003
  * Günther & Daylan (2021) —. 2021, ApJS, 254, 13, doi: 10.3847/1538-4365/abe70e
  * Hubbard (1984) Hubbard, W. B. 1984, Planetary interiors
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 435, 2152, doi: 10.1093/mnras/stt1435
  * Koch et al. (2010) Koch, D. G., Borucki, W. J., Basri, G., et al. 2010, ApJ, 713, L79, doi: 10.1088/2041-8205/713/2/L79
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161, doi: 10.1086/683602
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library.  http://ascl.net/1812.013
  * Lissauer (1995) Lissauer, J. J. 1995, Icarus, 114, 217, doi: 10.1006/icar.1995.1057
  * Lithwick et al. (2012) Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122, doi: 10.1088/0004-637X/761/2/122
  * Love (1911) Love, A. E. H. 1911, Some Problems of Geodynamics
  * Maciejewski et al. (2016) Maciejewski, G., Dimitrov, D., Fernández, M., et al. 2016, A&A, 588, L6, doi: 10.1051/0004-6361/201628312
  * Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, ApJ, 580, L171, doi: 10.1086/345520
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics
  * NASA Exoplanet Archive (2019) NASA Exoplanet Archive. 2019, Confirmed Planets Table, IPAC, doi: 10.26133/NEA1
  * Ni, Dongdong (2018) Ni, Dongdong. 2018, A&A, 613, A32, doi: 10.1051/0004-6361/201732183
  * Saillenfest et al. (2019) Saillenfest, M., Laskar, J., & Boué, G. 2019, A&A, 623, A4, doi: 10.1051/0004-6361/201834344
  * Seager & Hui (2002) Seager, S., & Hui, L. 2002, ApJ, 574, 1004, doi: 10.1086/340994
  * Storch & Lai (2013) Storch, N. I., & Lai, D. 2013, Monthly Notices of the Royal Astronomical Society, 438, 1526, doi: 10.1093/mnras/stt2292
  * Van Eylen & Albrecht (2015) Van Eylen, V., & Albrecht, S. 2015, ApJ, 808, 126, doi: 10.1088/0004-637X/808/2/126
  * Ward (1975) Ward, W. R. 1975, Science, 189, 377, doi: 10.1126/science.189.4200.377
  * Ward & Hamilton (2004) Ward, W. R., & Hamilton, D. P. 2004, AJ, 128, 2501, doi: 10.1086/424533
  * Winn et al. (2007) Winn, J. N., Holman, M. J., Henry, G. W., et al. 2007, AJ, 133, 1828, doi: 10.1086/512159
  * Zhu et al. (2014) Zhu, W., Huang, C. X., Zhou, G., & Lin, D. N. C. 2014, ApJ, 796, 67, doi: 10.1088/0004-637X/796/1/67

## appendix: The Effect Of Projection

When observing an oblate planet what is seen is a two-dimensional projection,
which is itself still an ellipse which can be characterised by a projected
oblateness and projected obliquity. We consider the amount by which the true
oblateness is altered if we assume a planet to have any uniform orientation in
space. Note that it is only possible for projection to decrease the oblateness
factor. When considering re-orientations of the planet, which has two long
major axes A and a short major axis $B=A(1-f)$, the radii of the cross section
of any rotation must be larger than B and smaller than A. Namely the projected
short axis b must obey $B_{\mathrm{proj}}\geq B$ and the project long axis
$A_{\mathrm{proj}}\leq A$. Recalling the definition of oblateness we find

$f_{\mathrm{proj}}=1-\frac{B_{\mathrm{proj}}}{A_{\mathrm{proj}}}\leq
1-\frac{B_{\mathrm{proj}}}{A}\leq 1-\frac{B}{A}\leq f$ (14)

When considering random rotations, in figure 10 we show the average
transformation of a uniform distribution of oblateness values. As expected we
see a noticeable shift towards shorter oblateness values. In this case we have
chosen the two orientation angles of the planet to be uniformly distributed
across the unit sphere. Analysing a large number of random projections, we
find that any given oblateness value is on average decreased by a factor of
about 20%.

Figure 10: The effect of projecting a three-dimensional ellipsoid onto the sky
plane, when its rotational axis is uniformly distributed across the unit
sphere. As expected there is general shift towards lower oblateness values.
This implies an additional level of difficulty when attempting to measure
oblateness, since even if a planet were to have a high oblateness it would
present itself as closer to spherical. Utilising a more informed model for the
orientation of planets would provide a more accurate description of the
transformation, however this would still be plagued by uncertainty about the
true underlying oblateness distribution.

We may also ask, given a measurement from say an MCMC fit to a planet in which
one retrieves a best fit value along with an uncertainty for the (projected)
oblateness and obliquity, what possible underlying distribution of the true
shape properties would match the observed values. For a choice of the
azimuthal angle, there is a unique mapping between the projected and un-
projected values, which can be obtained using equations 3a,3b,3c. The results
of this are shown in figure 11. For an azimuthal angle of 90 degrees we find
that the two distributions are the same. In this case the rotational axis of
the planet is tipped towards the observer, thus the two semi-major axes seen
are always the true long and short axis (implying
$f_{\mathrm{proj}}=f_{\mathrm{true}}$) and the obliquity angle is similarly
observed to be the same as the un-projected value.

Figure 11: The solid black cross and oval represent the 1$\sigma$ measurements
of a potential planet’s (projected) oblateness. The coloured patches represent
various distributions for the true (i.e. un-projected) values of oblateness
and obliquity which would match the projected observation, for different
choices of the azimuthal angle.