# Effects of tidal deformation on planetary phase curves

Babatunde Akinsanmi, Monika Lendl, Gwenael Boué, Susana C. C. Barros

► To cite this version:

Babatunde Akinsanmi, Monika Lendl, Gwenael Boué, Susana C. C. Barros. Effects of tidal deformation on planetary phase curves. Astronomy & Astrophysics - A&A, 2024, 682, 10.1051/0004-6361/202347739. insu-04822215

HAL Id: insu-04822215

https://insu.hal.science/insu-04822215v1

Submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effects of tidal deformation on planetary phase curves

Babatunde Akinsanmi<fnref n="1" />, Monika Lendl<fnref n="1" />, Gwenael Boué<fnref n="2" />, and Susana C. C. Barros<fnref n="3" />,<fnref n="4" />

Received 17 August 2023 / Accepted 4 October 2023

## ABSTRACT

With the continuous improvement in the precision of exoplanet observations, it has become feasible to probe for subtle effects that can enable a more comprehensive characterization of exoplanets. A notable example is the tidal deformation of ultra-hot Jupiters by their host stars, whose detection can provide valuable insights into the planetary interior structure. In this work we extend previous research on modeling deformation in transit light curves by proposing a straightforward approach to account for tidal deformation in phase curve observations. The planetary shape is modeled as a function of the second fluid Love number for radial deformation $h_{2f}$. For a planet in hydrostatic equilibrium, $h_{2f}$ provides constraints on the interior structure of the planet. We show that the effect of tidal deformation manifests across the full orbit of the planet as its projected area varies with phase, thereby allowing us to better probe the planet's shape in phase curves than in transits. Comparing the effects and detectability of deformation by different space-based instruments (JWST, HST, PLATO, CHEOPS, and TESS), we find that the effect of deformation is more prominent in infrared observations where the phase curve amplitude is the largest. A single JWST phase curve observation of a deformed planet, such as WASP-12 b, can allow up to a $17\sigma$ measurement of $h_{2f}$ compared to $4\sigma$ from transit-only observation. This high-precision $h_{2f}$ measurement can constrain the core mass of the planet to within 19% of the total mass, thus providing unprecedented constraints on the interior structure. Due to the lower phase curve amplitudes in the optical, the other instruments provide $\le 4\sigma$ precision on $h_{2f}$ depending on the number of phase curves observed. We also find that detecting deformation from infrared phase curves is less affected by uncertainty in limb darkening, unlike detection in transits. Finally, the assumption of sphericity when analyzing the phase curve of deformed planets can lead to biases in several system parameters (radius, dayside and nightside temperatures, and hotspot offset, among others), thereby significantly limiting their accurate characterization.

**Key words.** planets and satellites: interiors – planets and satellites: fundamental parameters – planets and satellites: atmospheres – techniques: photometric

## 1. Introduction

In ultra-hot-Jupiter systems, the tidal interaction between the star and the planet results in the deformation of both bodies into ellipsoidal shapes. The deformation of the star leads to a periodic modulation of the observed flux as the stellar tide, raised by the planet, rotates with the planet's orbital phase (Russell 1912; Morris 1985; Mislis et al. 2012; Esteves et al. 2013). This periodic signal in the light curve is referred to as stellar ellipsoidal variation, a phenomenon that has been empirically confirmed in several hot-Jupiter systems (see, e.g., Shporer et al. 2019; Wong et al. 2021; Parviainen et al. 2022). Similarly, the stellar tidal force deforms the planet along the axis pointing toward the star. The tidal deformation of the planet is challenging to measure during transit as the elongated axis is predominantly aligned with the line of sight, thereby resulting in a nearly circular sky projection of the planet. Nonetheless, slight deviations from the standard spherical planet transit light curve can be observed as the shape and cross-sectional area projected by the deformed planet varies with phase, thereby allowing the detection of deformation in high-precision light curves (Akinsanmi et al. 2019; Hellard et al. 2019; Barros et al. 2022). Since the projected area of the deformed planet varies throughout its orbit around the host star, it must also impact the observed orbital phase curve.

Although many effects are accounted for, including stellar ellipsoidal variation, when analyzing phase curve observations (see, e.g., Morris 1985; Loeb & Gaudí 2003; Esteves et al. 2013; Shporer 2017), the planet shape is usually assumed to be spherical due to the expected low amplitude of tidal deformation (Bourrier et al. 2020). However, with the precision of present and upcoming instruments, the subtle effects of deformation may now be detectable in phase curves. Since the deformation of a planet depends on its interior structure, its detection will offer valuable insight into the composition of the planet (Kellermann et al. 2018; Hellard et al. 2019; Akinsanmi et al. 2019). Furthermore, at high levels of precision, disregarding the effect of planetary deformation may bias interpretations of phase curves of ultra-hot Jupiters (UHJs).

In this paper, we drop the assumption of sphericity, but instead model the deformed planet shape as a triaxial ellipsoid to investigate the contribution of planet tidal deformation to phase curve observations of ultra-hot Jupiters. In Sect. 2 we model the shape of a deformed planet and describe a numerical and analytical approach to account for tidal deformation in phase curves of UHJs. In Sect. 3 we characterize the signature of tidal deformation across the full planetary orbit. We also compare the detectability of deformation in transit and phase curve observations for several space-based observatories. In Sect. 4 we illustrate how the measurement of tidal deformation can be used

<footnote n="1">¹ Observatoire astronomique de l'Université de Genève, chemin Pegasi 51, 1290 Versoix, Switzerland
e-mail: tunde.akinsanmi@unige.ch</footnote>

<footnote n="2">² IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote n="3">³ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="4">⁴ Departamento de Fisica e Astronomia, Faculdade de Ciencias, Universidade do Porto, Rua Campo Alegre, 4169-007 Porto, Portugal</footnote>
to constrain the core mass fraction of a planet, while in Sect. 5 we investigate biases in the parameters derived from phase curve analysis when sphericity is assumed for a deformed planet. We summarize our main conclusions in Sect. 6.

## 2. Modeling tidal deformation

### 2.1. Deformation and Love number

The deformation of a planet in response to stellar tidal forces depends on its interior structure and can be parameterized by the second-degree fluid Love number for radial displacement $h_{2f}$. In hydrostatic equilibrium, the giant planet deforms as a fluid body, and we have that $h_{2f} = 1 + k_{2f}$, where $k_{2f}$ is the second-degree Love number for potential that measures the distribution of mass within the planet<fnref n="1" /> (Love 1911; Kellermann et al. 2018). Since the Love number $h_{2f}$ depends on the distribution of mass within the planet, it measures the degree of central condensation of a body, thereby allowing us to constrain the interior structure (Kramm et al. 2011, 2012). The value of $h_{2f}$ ranges from 0 to 2.5, although certain nonequilibrium processes in the interior and nonlinear perturbative effects can result in values outside this range (Ragozzine & Wolf 2009; Wahl et al. 2021). A body with a fully homogeneous mass distribution (constant density) will have the maximum $h_{2f}$ value of 2.5, whereas a highly differentiated body with most of its mass condensed in a core will have a lower value ($h_{2f} \approx 1$ for main sequence stars; Ragozzine & Wolf 2009). More massive objects tend to have smaller $h_{2f}$ values since they are more compressible, and thus more centrally condensed (Leconte et al. 2011a). Measurements of $h_{2f}$ can therefore be used to infer the presence and mass of a planetary core. For instance, Saturn's lower $h_{2f}$ of 1.39 (Lainey et al. 2017) is indicative of a higher core mass fraction than Jupiter's $h_{2f} = 1.565$ (Durante et al. 2020). Similarly, limits on the Love number of HAT-P-13 b, obtained from its eccentricity due to its unique orbital configuration with a highly eccentric outer planet, allowed constraints to be placed on the maximum core mass of the planet (Batygin et al. 2009; Kramm et al. 2012; Buhler et al. 2016; Hardy et al. 2017). However, these works derive different limits on the Love number of the planet, highlighting the difficulty in accurately measuring the Love number of exoplanets. Therefore, measuring the shape of a tidally deformed giant planet from its light curve allows $h_{2f}$ to be estimated, which in turn provides direct constraints on the interior structure of the planet.

### 2.2. Shape model

The shape of a deformed planet is described here as a triaxial ellipsoid with semi-axes $[r_1, r_2, r_3]$ where $r_1$ is the planet radius oriented along the star-planet (substellar) axis, $r_2$ is along the orbital direction (dawn-dusk axis), and $r_3$ is along the polar axis. The volumetric radius of the ellipsoid is given by $R_v = (r_1 r_2 r_3)^{1/3}$. According to Correia (2014), the axes of the ellipsoid are related by $r_2 = R_v(1 - 2q/3)$, $r_1 = r_2(1 + 3q)$, and $r_3 = r_2(1 - q)$, where $q$ is an asymmetry parameter given by

$$q = \frac{h_{2f}}{2Q_M} \left(\frac{R_v}{a}\right)^3 \quad (1)$$

and depends on the Love number $h_{2f}$, the planet-to-star mass ratio $Q_M = M_p/M_\star$, and the orbital semimajor axis $a$.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Deformed planet phase curve components. (a) Phase-dependent projected area of an ellipsoid normalized to the projected area at phase zero (mid-transit). (b) Schematic of the projected ellipsoidal planet shape at different phases and the illuminated regions. (c) Atmospheric phase variation of a spherical (solid) and ellipsoidal (dashed) planet calculated analytically using Eq. (11). The cases with (green) and without (black) phase offset are shown. (d) Difference between the ellipsoidal and spherical planet atmospheric phase variations for the cases with and without phase offset. This represents the contribution of deformation to the planetary signal. For comparison, the stellar ellipsoidal variation $F_{EV}$ with a peak-to-peak amplitude of 100 ppm is also shown as the brown dotted curve. The cyan shaded regions in all panels represent the in-transit and in-eclipse phases.</figcaption></figure>

The nonspherical shape of the ellipsoidal planet implies that it projects a varying cross-sectional area as it rotates with orbital phase. The projected area as a function of orbital phase angle ($\phi = 2\pi \times \text{phase}$) is given (e.g., by Leconte et al. 2011a) as

$$A(\phi) = \pi \sqrt{r_1^2 r_2^2 \cos^2 \phi + r_3^2 \sin^2 \phi (r_1^2 \sin^2 \phi + r_2^2 \cos^2 \phi)}, \quad (2)$$

where $i$ is the orbital inclination of the planet. The projected area of the ellipsoid is plotted in Fig. 1a for a planet with the size, mass, and $h_{2f}$ of Jupiter orbiting a Sun-like star at a distance of 3 stellar radii. The area is normalized to the projected area at mid-transit. We see that the ellipsoid projects its minimum area at mid-transit and mid-eclipse (phases 0 and 0.5), while the maximum is at quadrature where the projected area is 9% larger than at mid-transit. As shown in previous studies (e.g., Correia 2014; Akinsanmi et al. 2019; Hellard et al. 2019), the varying projected shape and area of the ellipsoidal planet modifies the transit light curve. This allows the measurement of the planet's shape and its Love number, but requires high-precision transit data. Barros et al. (2022) were only able to detect the deformation of WASP-103 b by combining transit observations from several space-based instruments for improved precision. They obtained a $3\sigma$ significant measurement of $h_{2f}$. Hellard et al. (2020) obtained

<footnote n="1">¹ Calculations of the different Love numbers can be found in Sabadini & Vermeersen (2004) and Padovan et al. (2018).</footnote>
a less significant measurement for WASP-121 b using only a few transit observations from HST. Therefore, significantly detecting tidal deformation from transit observations remains quite challenging. We can further extend the previous works by probing for the varying area of the deformed planet in full-orbit phase curve observations, opening the potential for a more robust measurement of $h_{2f}$ that facilitates improved modeling of planetary interior structures.

## 2.3. Phase curve model

The phase curve model is a combination of planetary and stellar signals, which we summarize in this section with a modification introduced to account for planet deformation. The phase curve model is composed of transit $F_{tra}$, occultation $F_{occ}$, planet atmospheric phase variation $F_p$, and stellar phase variation $F_\star$ signals, and is given as a function of $\phi$ as

$$F(\phi) = F_{tra} \times F_{\star}(\phi) + F_{occ} \times F_{p}(\phi), \quad (3)$$

where the stellar phase variation is composed of the stellar ellipsoidal variation $F_{EV}$ (Morris 1985) and Doppler boosting $F_{DB}$ (Loeb & Gaudí 2003) signals as

$$F_{\star}(\phi) = 1 + F_{EV} + F_{DB}. \quad (4)$$

The components of Eqs. (3) and (4) are described in the following sections.

### 2.3.1. Transit and occultation

The planet's transit in front of the star and the occultation (secondary eclipse) behind the star can be modeled by several available tools. Here we use the publicly available ellc transit tool (Maxted 2018), where our transit model $F_{tra}$ adopts the power-2 limb-darkening law (Hestroffer 1997) with $c_{LD}$ and $\alpha_{LD}$ as limb-darkening coefficients (LDCs). The power-2 law has been shown to outperform other two-parameter laws in modeling intensity profiles generated by stellar atmosphere models (Morello et al. 2017; Claret & Southworth 2022). The occultation model, $F_{occ}$, is the same as the transit, but without limb darkening and normalized to a depth of unity.

To model the transit signal of a deformed planet, we followed the implementation of Akinsanmi et al. (2019) that incorporates the Correia (2014) ellipsoidal shape model (Sect. 2.2) into ellc. Therefore, our generated transit light curve takes into account the varying shape and cross-sectional area of the ellipsoid during transit. In addition to the usual transit parameters, the ellipsoidal transit model takes $h_{2f}$ and $Q_M$ as inputs, while replacing the typical spherical planet radius $R_p$ by the volumetric radius $R_v$ (Eq. (1)). Spherical planet transit and occultation signals can also be generated by setting $h_{2f} = 0$, which makes $R_v = R_p$.

### 2.3.2. Ellipsoidal variation

As previously mentioned, an ellipsoidal variation (EV) signal can be observed in phase curves of UHJs as a result of stellar deformation by a massive planet (Morris 1985; Mislis et al. 2012). Detailed expressions for the flux variations due to EV can be found in Esteves et al. (2013) and Csizmadia et al. (2023). However, the EV signal can be approximated by a cosine function at the first harmonic of the orbital period (Shporer 2017), so it has two peaks (at phases 0.25 and 0.75) during an orbit. It is given as

$$F_{EV} = A_{EV}(1 - \cos 2\phi), \quad (5)$$

where the semi-amplitude $A_{EV}$ depends on the planet-to-star mass ratio $Q_M$, the scaled semimajor axis $a/R_*$, and the inclination $i$ as

$$A_{EV} = \alpha_{EV} Q_M \left(\frac{R_*}{a}\right)^3 \sin^2 i, \quad (6)$$

where $\alpha_{EV}$ is a factor that depends on the linear limb-darkening coefficient and gravity-darkening coefficient.

### 2.3.3. Doppler beaming

This is an orbital photometric modulation caused by the motion of the star around the system's center of mass. The motion results in a Doppler shift of stellar light leading to a variation in the amount of photons observed in a particular passband (Loeb & Gaudí 2003; Shporer 2017). The signal is described by a sine curve at the orbital period, and so it has one peak (at phase 0.25) during an orbit. It is given by

$$F_{DB} = A_{DB} \sin(\phi), \quad (7)$$

where the semi-amplitude $A_{DB}$ is given by

$$A_{DB} = \alpha_{DB} 4 \frac{K}{c}. \quad (8)$$

Here $K$ is the RV semi-amplitude, $c$ is the speed of light, and $\alpha_{DB}$ is a coefficient of order unity that depends on the target's spectrum in the observed passband. For most planets with sub-km s⁻¹ RV amplitudes, $A_{DB}$ is very small (<5 ppm), and so Doppler beaming has little effect on the phase curve signal.

### 2.3.4. Planetary atmospheric phase variation

As a planet orbits its star, the planetary flux emanating from different fractions of the projected hemisphere contributes to the total observed flux from the system, resulting in a variation in the observed flux with phase. The planetary flux is composed of reflected light and thermal emission from the atmosphere, but these components are difficult to disentangle in observations as they produce variations that are similar in shape. However, since planets susceptible to tidal deformation are strongly irradiated, most of the observed flux is expected to be from thermal emission rather than reflected light. Here, we calculate the total planetary phase variation using numerical and analytical methods in order to investigate the contribution of planetary tidal deformation.

#### Numerical calculation

To numerically simulate the planetary phase variation signal, we adopted a pixelation approach where the planet (ellipsoidal or spherical) is projected on a two-dimensional Cartesian grid (XY). The projected shape of the ellipsoidal planet on the grid is an ellipse with semi-axes, $r_x$ and $r_y$, that varies with orbital phase and can be derived from Eq. (2). On the other hand, the projected shape of the spherical planet is a circle with radius $r_x=r_y$, which remains constant with orbital phase.

At each phase a flux map is assigned to the projected planet. Different flux maps can be adopted (see, e.g., Louden & Kreidberg 2018, and references therein). Here we define a flux map given as a function of co-latitude ($\theta$) and longitude ($\lambda$) of the visible hemisphere as

$$\mathcal{F}(\theta, \lambda) = A_0 + A_1 \sin \theta \cos(\lambda + \delta), \quad (9)$$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Schematic of the numerical calculation of planetary phase variation for an ellipsoidal and a spherical planet. **Top:** planet flux map for an ellipsoidal planet at transit (**left**), quadrature (**middle**), and eclipse (**right**) phases. The projected area of the planet changes with the number of pixels at each phase given as *npix*. **Middle:** same as top panels, but for a spherical planet. **Bottom:** planet flux variation obtained from integrating the flux map at the different phases. The black and red circles indicate the integrated flux value at each illustrated phase for the ellipsoidal and spherical planets, respectively.</figcaption></figure>

which simulates the simple sinusoidal phase function typically
used to model the total brightness variation of planets (e.g.,
Cowan et al. 2012; Kreidberg et al. 2018; Shporer et al. 2019;
Bell et al. 2019; Wong et al. 2021; Zieba et al. 2022). Here θ
varies in the range 0–180° between the poles, while λ varies from
λ_c − 90° to λ_c + 90° with λ_c denoting the central longitude of
the planet at a particular orbital phase (see Fig. 2). The longi-
tudes are defined such that the substellar point is on λ = 0° and
the antistellar point is on λ = 180°. The hotspot offset is denoted
by δ. The coefficients A₀ and A₁ define the flux value of pix-
els at the dawn–dusk limbs of the planet, and the flux difference

between these limbs and the center at eclipse, respectively. For
illustration, we set A₀ = 0 to simulate a total nightside flux of
zero, while we set A₁ = 1 so as to have a flux value of 1 at the
substellar point. The planet's phase variation Fp(φ) can then be
obtained by integrating the flux map across the projected surface
at each phase.

Assuming a planet with the size, mass, and h₂f of Jupiter
with a/R* = 3, Fig. 2 shows the flux map for the ellipsoidal
and spherical planet cases at transit, quadrature, and eclipse
phases. We note the larger projected area of the ellipsoidal
planet at quadrature where there are more brightly emitting
pixels compared to the spherical planet. The generated phase variation signals are also shown, where we see that the varying projected area of the ellipsoidal planet leads to additional flux between transit and eclipse since the emitting surface is larger at these phases. The maximum difference between the ellipsoidal and spherical planet occurs slightly after quadrature where the observer sees the largest illuminated or emitting surface of the ellipsoidal planet. Meanwhile, by construction, the spherical planet signal is equivalent to a cosine function. Numerically computing the planetary contribution to the phase curve can be computationally intensive, especially when fitting to observations. Therefore, we consider an analytical solution.

## Analytical calculation

Following previous phase curve studies (e.g., Cowan et al. 2012; Shporer et al. 2019; Wong et al. 2021), we analytically modeled the total atmospheric phase variation with a cosine function between the minimum flux, $F_{min}$, and maximum flux $F_{max}$, as

$$F_{p}(\phi) = (F_{max} - F_{min}) \frac{1 - \cos(\phi + \delta)}{2} + F_{min}. \quad (10)$$

The dayside flux $F_d$ (i.e., occultation depth) and nightside flux $F_n$ of the planet are derived as the value of $F_p(\phi)$ at $\phi = \pi$ and $\phi = 0$, respectively, while the semi-amplitude of the atmospheric variation $A_{atm}$ is $(F_{max} - F_{min})/2$.

As seen in the numerical simulation, the main effect introduced by tidal deformation in phase curves is geometric owing to the changing projected area of the ellipsoid. Therefore, we account for tidal deformation by multiplying the planet's phase variation (Eq. (10)) by the normalized phase-dependent projected area of the ellipsoid<fnref n="2" /> (Eq. (2)):

$$F_{p}(\phi)^{\text{def}} = F_{p}(\phi) \frac{A(\phi)}{A(0)}. \quad (11)$$

We note that other phase variation models (e.g., the kinematic temperature model of Zhang & Showman 2017, or other analytical functions, as given in Csizmadia et al. 2023) can be adopted for $F_p(\phi)$ and the same operation can be performed to incorporate tidal deformation. The analytical model produces an equivalent planetary phase variation signal to that obtained numerically (see Fig. A.1), thereby validating the simplified model that we use hereafter.

Figure 1c shows the atmospheric phase variations with zero phase offset for a spherical planet (green solid curve) and an ellipsoidal planet (green dashed curve) both with $F_d = 3600$ ppm<fnref n="3" />. The phase variation signals with eastward hotspot offset of 20° (black curves) are also shown. Figure 1d shows the difference between the phase variations of the ellipsoidal and spherical planets (i.e., $F_p(\phi)^{\text{def}} - F_p(\phi)$), highlighting the contribution of tidal deformation. In the case without hotspot offset (green curve), the contribution is symmetric about the mid-eclipse phase, peaking between quadrature and eclipse on both sides with an amplitude of ~200 ppm. In the case of a hotspot offset (black curve), the deformation contribution is asymmetric with some of its power shifted to the peak before

<figure><img src="image_4.png" /><figcaption>Fig. 3. Ellipsoidal planet phase curves with different values of $h_{2f}$, where $h_{2f} = 0$ is equivalent to a spherical planet. The inset zooms in on the phase curve region after transit.</figcaption></figure>

eclipse giving a higher amplitude of ~250 ppm before eclipse and 150 ppm after eclipse. The plot also shows the stellar ellipsoidal variation signal $F_{EV}$ with an arbitrarily chosen peak-to-peak amplitude of 100 ppm. The EV signal peaks at quadrature (i.e., spaced by 0.5 in phase) making it distinguishable from the planetary deformation signal, which has peaks at different phases with shorter spacing between them.

Our total phase curve implementation accounts for deformation by replacing $F_{tra}$, $F_{occ}$, and $F_p(\phi)$ in Eq. (3) by their ellipsoidal planet equivalents<fnref n="4" />. Figure 3 shows the ellipsoidal planet phase curve for different values of $h_{2f}$ with $h_{2f} = 0$ corresponding to a spherical planet. We see that the transit depth varies for different values of $h_{2f}$ since larger $h_{2f}$ implies more deformation, which causes the ellipsoidal planet to project a smaller cross section of its shape during transit. However, larger values of $h_{2f}$ increase the projected area of the ellipsoid toward quadrature, thereby causing an increase in the flux received from the planet between transit and occultation.

# 3. Signature and detectability of tidal deformation

Tidal deformation is only significant for UHJs orbiting close to the Roche limit of their stars, which limits the detectability to only a handful of confirmed planets. Previous studies (e.g., Correia 2014; Akinsanmi et al. 2019; Hellard et al. 2019; Berardo & de Wit 2022a) have identified WASP-12 b, WASP-103 b, and WASP-121 b as some of the most promising targets to measure deformation. For the analyses presented below, we adopt the parameters of WASP-12 b as representative of such UHJs.

## 3.1. The signature of deformation

The detectability of deformation requires that the observed light curve of a deformed planet is distinct such that it cannot be accurately reproduced by a spherical planet model. From visual inspection of Fig. 3, one can imagine that most of the differences between the spherical and deformed planet phase curves can perhaps be explained by a modification of some of the spherical

<footnote n="2">² This is done only for the out-of-transit and occultation phases since the projected shape and area of the deformed planet is already accounted for in the transit and occultation models (Sect. 2.3.1).</footnote>

<footnote n="3">³ Similar to *Spitzer* 3.6 µm and 4.5 µm occultation depth measurements for a planet with a dayside temperature of ~2700 K orbiting a Sun-like star.</footnote>

<footnote n="4">⁴ The code is publicly available at https://github.com/tundeakins/tidef</footnote>
<table><caption>**Table 1.** Parameters adopted in our optical and NIR phase curve simulations of WASP-12 b.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Optical/TESS<sup>(†)</sup></th><th>NIR/PRISM<sup>(‡)</sup></th></tr></thead><tbody><tr><td>R<sub>p</sub> (R<sub>*</sub>)</td><td>Scaled spherical radius</td><td>0.1160</td><td>0.1168</td></tr><tr><td>R<sub>v</sub> (R<sub>*</sub>)</td><td>Scaled volumetric radius<sup>(§)</sup></td><td>0.1222</td><td>0.1232</td></tr><tr><td>a/R<sub>*</sub></td><td>Scaled semimajor axis</td><td>3.061</td><td>3.061</td></tr><tr><td>b</td><td>Impact parameter</td><td>0.344</td><td>0.344</td></tr><tr><td>c<sub>LD</sub></td><td rowspan="2">Power-2 LDCs<sup>(§)</sup></td><td>0.6371</td><td>0.2862</td></tr><tr><td>α<sub>LD</sub></td><td>0.6175</td><td>0.4183</td></tr><tr><td>F<sub>d</sub>/F<sub>*</sub> (ppm)</td><td>Dayside flux</td><td>466</td><td>3306</td></tr><tr><td>F<sub>n</sub>/F<sub>*</sub> (ppm)</td><td>Nightside flux</td><td>0</td><td>500</td></tr><tr><td>δ (°)</td><td>Phase offset</td><td>6.2</td><td>20</td></tr><tr><td>A<sub>EV</sub> (ppm)</td><td>EV semi-amplitude</td><td>65</td><td>60</td></tr><tr><td>A<sub>DB</sub> (ppm)</td><td>DB semi-amplitude</td><td>2.34</td><td>1.18</td></tr><tr><td>h<sub>2f</sub></td><td>Love number<sup>(§)</sup></td><td></td><td>1.565</td></tr><tr><td>Q<sub>M</sub></td><td>Planet-to-star mass ratio<sup>(§)</sup></td><td></td><td>0.00098</td></tr></tbody></table>

**Notes.** (†) The TESS parameter values are the derived values from the analysis of four sectors of WASP-12 TESS observations Wong et al. (2022). (§) These values were obtained using LDTk (Parviainen 2018) in the different passbands. (‡) Values for the PRISM parameters R<sub>p</sub> and F<sub>d</sub>/F<sub>*</sub> are obtained by integrating the modeled transmission (Stevenson et al. 2014) and emission spectra (Hooton et al. 2019) across the PRISM wavelength range, A<sub>EV</sub> and A<sub>DB</sub> are estimated using Eqs. (6) and (8), respectively, while δ is the derived value from the 2013 *Spitzer* 4.5 µm phase curve of WASP-12b (Bell et al. 2019). (§) These parameters are specific to the ellipsoidal planet model where we assume the same h<sub>2f</sub> as Jupiter, R<sub>v</sub> is the volumetric radius that produces the same transit depth as the spherical planet of radius R<sub>p</sub>, and Q<sub>M</sub> is taken from Collins et al. (2017).

planet model parameters (e.g., radius, hotspot offset). However, in actual observations, the value of these parameters would be initially unknown, but determined by performing a fit to the data. Therefore, the detectable signature of deformation is the residuals obtained from fitting a deformed planet observation (transit or phase curve) with a spherical planet model. The amplitude and timescale of the residual anomalies indicate the prominence of deformation in the data and hence its detectability.

To estimate and compare the signature of deformation in transit and phase curve observations, we assume broadband observations of WASP-12 b in the optical (e.g., with TESS, CHEOPS, or PLATO) and band-integrated observations in the near-infrared (NIR; e.g., with JWST or HST). We compared the deformation signatures by simulating ideal transit and phase curve observations (but no uncertainties) of a deformed planet and performing a least-squares fit<fnref n="5" /> using a spherical planet model. We selected TESS for the representative simulation of the optical observations based on the available measurements and JWST/NIRSpec in PRISM mode for the NIR. Although we used TESS and NIRSpec-PRISM as examples of observing instruments, the approach is applicable to other instruments observing in the optical and NIR. The adopted parameters of the simulation are given in Table 1; they are all allowed to vary freely except for the LDCs and A<sub>DB</sub>, for which we used Gaussian priors. In the transit case, the out-of-transit baseline is kept short to minimize the effect of orbital phase variation on the fit, and we put strong priors on the parameters (F<sub>d</sub>, F<sub>n</sub>, δ, A<sub>EV</sub>) since they cannot be

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Fitting the phase curve and transit of a deformed planet with a spherical planet model. Left: simulated deformed planet phase curve of WASP-12 b in the TESS (blue) and NIRSpec-PRISM (red) bands. The spherical planet model fit is overplotted for each light curve. The bottom panel shows the residuals of the fit for each light curve. Right: same as the left panel, but for observation of the transit only. The residuals are shown in the bottom panel for each light curve with the root sum of squares (RSS) quoted.</figcaption></figure>

accurately estimated from the transit fit. For real observations, the phase variation in transit-only observations can be accounted for either by applying priors on these parameters from previous phase curve analyses and observations of the target or by detrending using a time-dependent function.

The results of the fits are shown in Fig. 4. The residuals of the transit fits display anomalies that represent the signature of deformation in transit light curves as described in Correia (2014) and Akinsanmi et al. (2019): transit ingress and egress anomalies due to the longer transit of the ellipsoidal planet, and in-transit anomalies due to varying projected shape and eclipsed area by the rotating ellipsoidal planet. The residuals of the phase curve fits show an additional anomaly across the orbit also due to the varying projected area of the ellipsoidal planet, which contributes additional flux at different phases that cannot be easily accounted for by the spherical planet model. There are also ingress and egress anomalies in the occultation due to the longer occultation duration of the ellipsoidal planet.

To quantify the deviation of the spherical model fits from the deformed planet simulated observations, we computed the root sum of squares (RSS) of the residuals. We find a greater deviation in the PRISM transit residuals than for TESS indicating a larger deformation signal in the NIR than the optical due to the reduced effect of limb darkening in the NIR that amplifies the ingress and egress anomalies. The phase curve residuals also exhibit larger deviations than the corresponding transit residuals (>2× for PRISM) due to the wider phase range of the phase curve deformation signature. This implies that the effect of deformation is more prominent in phase curves than in transit-only observations. Comparing the phase curve deformation signature between the instruments, we see that the out-of-transit amplitude is just as high as the in-transit amplitude for PRISM, whereas it is much smaller for TESS. The RSS of the TESS phase curve signature does not provide much improvement compared to the RSS of the transit due to the small contribution from the out-of-transit signature.

<footnote n="5">⁵ With LMFIT (https://lmfit.github.io/lmfit-py/; Newville et al. 2020).</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Detectability of deformation in the optical (TESS; left) and NIR (JWST/NIRSpec-PRISM; right). The horizontal dashed line indicates the simulated value of $h_{2f} = 1.565$. The red points show the retrieved value of $h_{2f}$ from phase curves and the $1\sigma$ uncertainty while the blue points show the retrieved $h_{2f}$ from transits. The gray dotted error bars show the $3\sigma$ uncertainty of the measurements. A similar plot for HST/WFC3 is shown in Fig. A.2.</figcaption></figure>

In summary, phase curves provide a longer phase range than transits, which allows the shape of the planet to be probed better as it rotates with phase. The longer phase range can also be more easily sampled to attain better precision, which facilitates the detection of the anomaly. We quantify the detectability of deformation from transits and phase curves in the next section.

## 3.2. Detectability of tidal deformation in transit and phase curve observations

Detecting tidal deformation implies measuring a statistically significant value of $h_{2f}$ from fitting transit and/or phase curve observations with the ellipsoidal planet model described in Sect. 2. We investigated the detectability of tidal deformation at optical and NIR passbands by performing injection-retrieval simulations. We simulated the transit and phase curve of deformed WASP-12 b at a cadence of two minutes using the parameters listed in Table 1. For ease of fitting, we did not include any phase variation in the simulated transit signals, assuming that this will be accounted for as previously mentioned.

In order to investigate the photometric precision required to detect deformation and assess the detectability by different instruments, we added random Gaussian noise of different levels to the simulated transit and phase curve, which we then fit with ellipsoidal and spherical planet models. The *dynesty* Python package (Speagle 2020) was used to sample the parameter space and derive parameter posteriors. All parameters for each model were allowed to vary freely except for the power-2 LDCs, for which we used Gaussian priors with arbitrary uncertainties of 0.01, $Q_M$ whose prior was based on derived values in Collins et al. (2017), and $A_{DB}$ whose prior was derived from Eq. (8).

For each noise level we generated 20 simulated observations by adding random realizations of the noise. We then performed the model fit to each observation. The posterior distributions from the fits with the different noise realizations were then merged. This method marginalizes over several possible noise realizations, and thus mitigates the impact of obtaining biased results based on a single randomly generated favorable or unfavorable noise realization.

The detectability of tidal deformation in the optical and NIR is summarized in Fig. 5. It shows the retrieved value of $h_{2f}$ from the fits to the simulated optical and NIR observations at different noise levels. For the optical simulations, the transit fits allow $h_{2f}$ to be recovered with $3\sigma$ significance for noise levels up to 70 ppm/2 min, while the phase curve fits extend this to a higher noise level of 85 ppm/2 min. At noise levels $\le 50$ ppm/2 min the phase curve fits provide up to 15% improvement in the precision of retrieved $h_{2f}$ compared to the transit fits. However, the improvement decreases below 5% at higher noise levels since the relatively small amplitude of the deformation signature outside of transit phases becomes indiscernible from the noise. For the NIR simulations, $h_{2f}$ is recovered at $3\sigma$ in the transit fits for noise levels up to 120 ppm/2 min. Across all added noise levels, the phase curve fits provide a significant improvement (76% at 50 ppm and $>30\%$ at 150 ppm) in the precision of retrieved $h_{2f}$ compared to the transit fits due to the long phase range and large amplitude deformation signature outside transit.

Our injection and retrieval simulations thus imply that detecting tidal deformation is more favorable at NIR passbands compared to the optical since we are able to measure $h_{2f}$ significantly at higher noise levels in the NIR, particularly for phase curve observations. As phase curves require more observing time than transits, we compared $h_{2f}$ retrieval from phase curves to the same observing time in transit observations. This allowed us to determine the most effective observing strategy for the precise measurement of $h_{2f}$. For the same observing time, we find that phase curves provide improved $h_{2f}$ precision over the combined transit by up to 10% in the optical and up to 50% for JWST.

### 3.2.1. Impact of limb darkening on detectability

The effect of limb darkening is most significant at ingress and egress, where the signature of tidal deformation also manifests strongly. Limb darkening is therefore capable of compensating for part of the deformation-induced signal, thus reducing detectability, particularly for transit observations (Hellard et al. 2019; Akinsanmi et al. 2019). Therefore, detecting deformation requires proper treatment of limb darkening.
We find in our injection-retrieval simulations that increas-
ing the width of our limb-darkening priors results in reduced
detection significance. For instance, when the widths of the
LDC priors are slightly increased from 0.01 to 0.015 in the
optical simulations, the noise level required for a 3σ detection
reduces from 70 to 50 ppm/2 min in the transit fits and from 85 to
55 ppm/2 min in the phase curve fits. The reduction in both cases
occurs because most of the detection power in the optical comes
from the transit phases where limb darkening acts. However, the
increased LDC prior width in the NIR causes only a slight reduc-
tion in noise level from 120 to 110 ppm/2 min for a 3σ transit
detection due to the weaker limb-darkening effect in the infrared.
The NIR phase curve fits remain largely unchanged, even with a
further increase in the LDC prior width, since the deformation
signature spans a wide phase range with a large out-of-transit
amplitude that allows a precise constraint of the deformation in
spite of limb darkening.

### 3.2.2. Disentangling deformation from atmospheric phenomena

It is possible for other atmospheric effects to confound the defor-
mation signal or make its detection challenging. For example,
highly irradiated planets may exhibit strong day-to-night temper-
ature gradients that can result in different chemical compositions
and scale heights between the two hemispheres (see, e.g., Pluriel
et al. 2020; Helling et al. 2020). This leads to limb asymme-
tries (Espinoza & Jones 2021) that can affect the detectability of
deformation in transit.

To investigate whether the effects of day-to-night inhomo-
geneities can mimic tidal deformation, we simulated the transit
of a planet with limb asymmetries using catwoman (Jones &
Espinoza 2020) and fit it with a deformed planet model. We
find that tidal deformation is unable to account for limb asymme-
tries since tidal deformation does not cause asymmetric transits.
Conversely, we also found that catwoman is unable to accurately
fit a deformed planet transit indicating that the two effects are
distinguishable.

Furthermore, the effect of dayside–nightside inhomo-
geneities is chromatic and is dependent on the compositions of
the two regions. In contrast, the Love number is an intrinsic
property of the planet that remains constant across wavelengths.
Therefore, we can differentiate between them by ensuring that
consistent $h_{2f}$ is derived across different spectral bins or pass-
bands. Therefore, spectroscopic observations with HST and
JWST will be very useful in this regard.

### 3.2.3. Detectability of deformation by different instruments

Figure 5 can be used to estimate the detectability of deformation
by different instruments observing phase curves or transits of
close-in giant planets in the optical and NIR. Here, we estimate
the detectability for observations with JWST, PLATO, TESS,
CHEOPS, and HST.

**JWST.** Using the JWST Exposure Time Calculator,<fnref n="6" /> we estimated a precision of 40–50 ppm/2 min for observations of WASP-12 b with different JWST instruments, for example NIRSpec-PRISM. At the similar noise level of 50 ppm/2 min, our NIR simulations recovered a highly significant detection of $h_{2f} = 1.52^{+0.094}_{-0.087}$ (17σ) from the phase curve fit and $h_{2f} = 1.55^{+0.32}_{-0.37}$ (4.2σ) from the transit fit. Figure 6 shows the ellipsoidal and

spherical planet model fits to the simulated PRISM observa-
tions with one realization of the 50 ppm noise. The residual of
the spherical planet fit shows a clear modulation induced by
tidal deformation that cannot be accounted for by the spher-
ical planet. Figure 6 further shows the contribution of planet
deformation (cf. Fig. 1d) to the total observed phase curve in
comparison to the best-fit stellar ellipsoidal variation signals
$F_{EV}$. We observe that the amplitude of $F_{EV}$ in the spherical
planet model fit increases (from the simulated value of 60 ppm
to 210 ppm) in an effort to account for the deformation signal of
the planet. However, it is unable to do so since the deformation
contribution peaks at different phases with a shorter inter-peak
interval compared to $F_{EV}$. This is still the case when a tidal phase
lag parameter is added to $F_{EV}$. This implies that the assump-
tion of sphericity for a deformed planet can lead to biases in
the estimation of other system parameters. In Sect. 5 we fur-
ther investigate the bias on different system parameters when
assuming sphericity for a deformed planet.

**PLATO.** For stars brighter than $V \approx 11$, PLATO will be capable of achieving a photometric precision of 34 ppm h⁻¹ (Rauer et al. 2014; Rauer & Heras 2018) equivalent to 186ppm/2 min. The observation strategy of PLATO consists of a long-pointing phase where two different fields are observed for two to three years each, and a short observation phase (known as step and stare), where each stare would last two to three months (Rauer & Heras 2018). Assuming a conservative scenario where WASP-12 is observed in the short observation mode for two months, this will amount to ~50 orbits of the planet leading to a precision of ~30 ppm/2 min. At this noise level in our optical simulations, we recovered a 4.3σ detection with $h_{2f} = 1.59^{+0.40}_{-0.37}$ from the phase curve fit and a 3.6σ detection with $h_{2f} = 1.54 \pm 0.42$ from the transit fit. Figure A.4 shows the best-fit models to the simulated PLATO observations. PLATO long-pointing observations of 2–3 yr of UHJs susceptible to deformation will provide exquisite photometric precision that will only be limited by stellar noise and limb-darkening modeling.

**TESS.** We estimated the precision of TESS from its actual observations of WASP-12 b for which it obtained a precision of ~1800 ppm/2 min (Wong et al. 2022). A total of 78 orbits of the planet were observed (across four TESS sectors), which improved the precision to ~200 ppm/2 min. At this noise level in our optical simulations, we obtained a ~2.3σ detection of $h_{2f}$ from the phase curve fit and a slightly lower ~2σ detection from the transit fit. Similar to Fig. 6, Fig. A.3 shows the best-fit models to the simulated TESS observations with 200 ppm/2 min noise. In order to detect deformation with higher significance, more TESS observations will be required to lower the noise level.

**CHEOPS.** Using the CHEOPS Exposure Time Calculator,<fnref n="7" /> we estimated a photometric precision of 450 ppm/2 min<fnref n="8" /> for observations of WASP-12 b. For a month-long observation of WASP-12 b (or 28 accumulated phase curves), a precision of 85 ppm/2 min can be attained. At this noise level in our optical simulations, we recovered a 3σ detection of $h_{2f}$ from the phase curve fit and 2.5σ from the transit fit. Figure A.5 shows the best-fit models to the simulated CHEOPS observations.

**HST.** Observations of WASP-12 b have been performed with different HST instruments. Using the G141 grism of HST/WFC3

<footnote n="6">⁶ https://jwst.etc.stsci.edu/</footnote>

<footnote n="7">⁷ https://cheops.unige.ch/pht2/exposure-time-calculator/</footnote>

<footnote n="8">⁸ Consistent with actual CHEOPS observation of WASP-12 b (Akinsanmi et al., in prep.).</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Spherical planet (SP) and ellipsoidal planet (EP) model fits to the simulated deformed planet signals with a PRISM noise level of 50 ppm/2 min. The top panel shows the simulated phase curve (**left**) and transit (**right**) observations with the fitted models overplotted. The phase curve fit also shows the contribution of deformation to the light curve (blue) and the $F_{EV}$ component of the EP and SP model fits. The middle panel shows the residuals from fitting a spherical planet model (cyan) and an ellipsoidal planet model to the simulations. **bottom panels** show the residuals from fitting a spherical planet model (red) and an ellipsoidal planet model to the simulations.</figcaption></figure>

Kreidberg et al. (2015) obtained a precision equivalent to ~110 ppm/2 min on three transit visits. Phase curve and transit simulations at this noise level result in a $4.3\sigma$ and $2.6\sigma$ detection of $h_{2f}$, respectively. We see that a noise level of 110 ppm for WFC3 constrains $h_{2f}$ better than 85 ppm from CHEOPS due to the larger NIR phase curve amplitude (Fig. A.6). Doubling the number of visits to six improves the detection significance in the phase curve and transit to $5.6\sigma$ and $3\sigma$, respectively.

We note that gaps within the data (particularly for HST and CHEOPS) will reduce the deformation detection significance estimated for the different instruments or increase the number of visits required to reach a certain detection significance. For example, even after combining three HST/WFC3-G141 transit visits of WASP-12 b in Kreidberg et al. (2015), the ingress and egress phases still have poor coverage, making it challenging to measure tidal deformation from these datasets within a reasonable number of HST visits. However, observations can be scheduled so that the combined visits attain better phase coverage. Instrumental and/or astrophysical noise (e.g., spacecraft jitter or stellar activity) can also impact the detection of deformation from the light curves. These will have to be mitigated while analyzing the data using methods such as polynomial decorrelations, Gaussian processes (Barros et al. 2012), or wavelets (Csizmadia et al. 2023).

## 4. Constraining the core mass from $h_{2f}$

As mentioned previously, measuring the Love number of a planet sheds insight into its interior structure. For instance, Baumeister et al. (2020) found that adding the Love number as input to interior structure models significantly decreases the degeneracy of the possible interior structure. Furthermore, simulations of the interior structure of low-mass planets by Baumeister & Tosi (2023) found that a Love number precision of $10\sigma$ can constrain the core and mantle size of an Earth analog to ~13% of its true

value. Here we show how the retrieved value of $h_{2f}$ from the transit and phase curve fits of the different instruments can be used to place constraints on the core mass fraction $M_c/M$ of the planet. To do this, we calculate the theoretical value of $h_{2f}$ as a function of $M_c/M$ using Eqs. (1) and (2) of Buhler et al. (2016) and adopting the radial density profile in Eq. (7) of Helled et al. (2014), which consists of a constant density core and a stratified density envelope.

Figure 7a shows the dependence of $h_{2f}$ on $M_c/M$. As $M_c/M$ increases the planet is more centrally condensed and less homogeneous, which causes $h_{2f}$ to decrease. We overplot the retrieved values of $h_{2f}$ from the transit and phase curve fits of the JWSTNIRSpec-PRISM simulations (Sect. 3.2.3). We see from the $1\sigma$ lower limits of $h_{2f}$ that the phase curve fit constrains the maximum core mass fraction of the planet as 0.19, whereas the larger $h_{2f}$ uncertainties from the transit fit give a maximum $M_c/M$ of 0.56. Figure 7b compares the maximum core mass fraction derived from the retrieved $h_{2f}$ from the phase curve and transit fit of the different instruments. We see that the JWST phase curve fit provides the tightest constraint on the core mass of the planet followed by the HST and PLATO phase curve fits, and then the JWST and PLATO transit fits. We see that the $<3\sigma$ $h_{2f}$ measurements from CHEOPS and TESS provides no constraint on the core mass with a maximum core mass being equal to the planetary mass. We note that although the measurements from the different instruments are still consistent with a zero core mass, more precise $h_{2f}$ measurements (e.g., from more than one JWST phase curve or PLATO long-pointing) will also allow the lower limit of $M_c/M$ to be determined, definitively constraining the planetary core and improving interior structure models.

## 5. Bias on derived system parameters

Parameters of the model adopted in fitting observations are often correlated with one another, making it possible for other
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Constraints on the core mass fraction of WASP-12 b derived from the retrieved *h*2f. **Left:** Constraints from JWST phase curve (red) and transit (blue). The black curve shows the theoretical Love number as a function of core mass fraction for a two-layer model using the formulations of Buhler et al. (2016) and Helled et al. (2014). **Right:** Upper limits on core mass fraction based on the precision of retrieved *h*2f from different instruments. The number of phase curve and transits vists used are JWST -1 visit, HST -3 visits, PLATO -50 visits, CHEOPS -28 visits, TESS -78 visits.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Parameter biases from fitting a spherical planet model to simulated phase curve of deformed planets of different *h*2f. The parameter names are given in Table 1. For each parameter, we compare the deviation of the fitted spherical planet parameters (fit) from the true input parameters of the deformed planet (truth) and scale this by the standard deviation of the fit (σfit). The colors show the biases for the different instruments.</figcaption></figure>

phase curve parameters to dampen or entirely mask the deformation signal depending on the noise level of the light curve. Assuming sphericity when fitting the phase curve of a deformed planet could introduce biases in the estimation of various other parameters. The deformation then manifests itself as an astrophysical source of systematic error in the measurement of the parameters.

We investigate the parameter biases by fitting a spherical planet model to the simulated TESS, PLATO, CHEOPS, HST, and PRISM deformed planet phase curves with increasing values of *h*2f up to 1.8. The noise level of each instrument estimated in Sect. 3.2.3 was added to the simulated light curves. For each

parameter we compare the obtained posterior to the true simulated value by calculating the number of standard deviations by which the median posterior value differs from the true value: (fit – truth)/σfit. This normalization allows us to easily compare the parameter deviation for each instrument.

The result of the fits is shown in Fig. 8. As expected, for *h*2f = 0 (spherical planet) the fitted parameters are close to (\<1σ) the true simulated values, leading to deviations close to zero. A zero deviation is also obtained for cases where the deformed planet phase curve is indistinguishable from that of a spherical planet, due to high noise level in the light curve compared to the deformation signal.
For all instruments we observe significant biases in the esti-
mation of R<sub>p</sub>/R<sub>⋆</sub>, a/R<sub>⋆</sub>, and b, which all show clear trends
with h<sub>2f</sub>. Since a tidally deformed planet projects only a small
cross section of its shape during transit, the best-fit spherical
planet model has a smaller R<sub>p</sub> than the volumetric radius R<sub>v</sub>
of the ellipsoid. Therefore, as h<sub>2f</sub> increases and the planet is
more deformed, it projects an even smaller cross-sectional area,
which requires a smaller spherical R<sub>p</sub> to match the transit depth.
We find that the derived R<sub>p</sub>/R<sub>⋆</sub> increasingly deviates from (i.e.,
underestimates) the simulated R<sub>v</sub>/R<sub>⋆</sub> as a function of h<sub>2f</sub>. At
the highest simulated h<sub>2f</sub> value of 1.8, R<sub>p</sub>/R<sub>⋆</sub> deviates by up
to 25σ for the TESS fit; 50σ for the CHEOPS fit; and >80σ
for the HST, PLATO, and PRISM phase curve fits, which have
the most precise radius measurements. When combined with
the mass estimate from radial velocity (RV) measurements, the
radius deviation can thus lead to an overestimation of the density
of the planet by more than 10% if sphericity is assumed (Burton
et al. 2014; Correia 2014; Barros et al. 2022). As the precision of
RV mass estimates increases, the bias in derived densities due to
tidal deformation will constitute a bottleneck in constraining the
composition of exoplanets susceptible to deformation (Berardo
& De Wit 2022b).

Both *a/R<sub>⋆</sub>* and *b* deviate at *h<sub>2f</sub>*=1.8 by >2.5σ for TESS, CHEOPS, and HST, while the higher precision phase curves of PLATO and PRISM deviate by >12σ. These parameters deviate from the true inputs in order to compensate for the longer transit duration of the deformed planet (Arkhypov et al. 2021). The limb-darkening parameters *c<sub>LD</sub>* and *α<sub>LD</sub>* also show biases as these parameters attempt to absorb the deformation signal concentrated around ingress and egress.

Biases are also observed in the dayside and nightside fluxes, F<sub>d</sub> and F<sub>n</sub>, which show opposing trends as a function of h<sub>2f</sub>. Both parameters act together to increase the amplitude A<sub>atm</sub> of the atmospheric phase variation in order to account for the additional flux from the deformed planet at phases with larger projected areas compared to the sphere (see Fig. 6). Biases on these parameters affect the temperatures derived for the hemispheres of the planet. The temperature of the dayside is overestimated and that of the nightside is underestimated, which results in a larger day-to-night contrast. The deviations in F<sub>d</sub> and F<sub>n</sub> are up to 20σ at h<sub>2f</sub>=1.8 for the PRISM fit where the amplitudes of these parameters are much higher than the light curve scatter. Converting the dayside and nightside fluxes to temperatures<sup>9</sup> (T<sub>d</sub> and T<sub>n</sub>), we find that their deviations lead to a 5.4σ overestimation of T<sub>d</sub> and a 13σ underestimation of T<sub>n</sub> for PRISM, and ≤3σ deviations for HST. For the optical phase curves with lower simulated values of F<sub>d</sub>=466 ppm and F<sub>n</sub>=0 ppm, deviations in the derived F<sub>d</sub> decrease to ~8σ (T<sub>d</sub> deviation of 2.2σ) for PLATO and <2σ for TESS and CHEOPS (T<sub>d</sub> deviations of <1σ), whereas deviations in F<sub>n</sub> are close to zero.

The hotspot offset δ also exhibits a positive deviation of up to 11σ at h<sub>2f</sub>=1.8 for the PRISM phase curves with the largest simulated offset of 20°. The larger eastward hotspot of the spherical planet attempts to compensate for the deformation-induced flux before eclipse. However, the optical simulations with small offsets of 6.2° showed only <1.5σ deviations. As mentioned in Sect. 3.2.3, the amplitude of ellipsoidal variation amplitude increases in an attempt to compensate for the deformation contribution. This leads to a positive deviation from the true value by up to 50σ at h<sub>2f</sub>=1.8 for the PRISM fit, ~10σ for the HST and PLATO fits, 5σ for the CHEOPS fit, and 2σ for the TESS

fit. There are no significant deviations in the derived transit time
and Doppler beaming amplitude.

In general, the PRISM phase curve shows the largest param-
eter deviations due to the high precision of JWST data and
the large amplitude deformation signal that the spherical planet
model parameters try to compensate for. Conversely, the TESS
phase curve shows the smallest parameter deviations due to the
high noise level of the simulated observations, which makes
the phase curve similar to that of a spherical planet and makes
detecting deformation challenging. Nevertheless, we see that
even for such low or nondetection cases, deformation can still
lead to significant deviations (ranging from 3 to 25σ) in some
of the derived parameters. Therefore, we recommend accounting
for planetary deformation when fitting the phase curves of UHJs
by multiplying the planet's atmospheric phase variation by the
phase-dependent projected area of the ellipsoid (Eq. (2)). This
introduces a single new parameter h2f to fit, if the mass ratio
is kept fixed to RV values, which allows marginalizing over the
possible shapes of the planet. In analyzing the HST and Spitzer
phase curves of WASP-103 b, Kreidberg (2018) accounted for the
planet's deformation by similarly including the phase-dependent
projected ellipsoid area in their model. However, instead of
fitting for shape parameters of the ellipsoid, the planet shape
was fixed using tabulated predictions in Leconte et al. (2011a)
assuming a planetary radius and age.

6. Conclusions

In this work we investigated the effects of tidal deformation on
the phase curves of UHJs orbiting close to the Roche limit of
their stars. We expanded on previous works (e.g., Leconte et al.
2011a; Burton et al. 2014; Correia 2014; Akinsanmi et al. 2019;
Hellard et al. 2019) studying tidal deformation in transit obser-
vations by introducing a simple approach to additionally account
for tidal deformation when modeling phase curve observations
of UHJs. Similar to Correia (2014) and Akinsanmi et al. (2019),
our model describes the planet as a triaxial ellipsoid parameter-
ized by the Love number h2f, thereby allowing us to estimate its
value from fitting observations. In addition to the signature of
deformation in transit caused by the varying projected shape and
area of the planet, we showed, using numerical and analytical
calculations, that the planet's deformation also manifests itself
outside of transit for the same reason. Therefore, the extended
phase coverage of the deformation signal in phase curves allows
the planet's shape to be probed better.

Using injection-retrieval simulations, we find that detecting
tidal deformation is more favorable at NIR passbands compared
to the optical as we are able to measure h<sub>2f</sub> significantly at higher
noise levels in the NIR. Furthermore, NIR phase curve observa-
tions significantly improve the deformation detection compared
to transit-only observations due to larger planetary phase curve
amplitudes in the NIR. For the precision expected from several
JWST instruments and modes, we find that the measurement of
h<sub>2f</sub> for WASP-12 b improves from a precision of ~4σ from one
transit-only analysis to 17σ from a single phase curve. Combin-
ing several transits equivalent to the phase curve time investment
still does not attain the h<sub>2f</sub> measurement significance of a phase
curve.

The high h<sub>2f</sub> precision obtainable from JWST phase curves
of relevant targets will provide unprecedented constraints on the
planet's interior structure (e.g., core mass fraction Buhler et al.
2016; Kramm et al. 2012). Additionally, spectrophotometric light
curves from HST and JWST will allow us to verify that the

<footnote n="9">⁹ Using a BT-Settl stellar model (Allard et al. 2012) for WASP-12 and a blackbody for the planet.</footnote>
derived $h_{2f}$ is not wavelength dependent, effectively disentangling tidal deformation from atmospheric phenomena that are chromatic.

Our results also show that detecting deformation from phase curves is relatively unaffected by the uncertainty in the limb-lightening profile as long as there is a high-amplitude deformation signature outside the transit phases, as obtainable for NIR observations. Finally, we showed that the assumption of sphericity when analyzing the phase curve of UHJs can lead to biases in several derived parameters, thereby significantly limiting the ability to properly characterize these targets (e.g., their radii, densities, dayside and nightside temperatures). We thus recommend accounting for deformation by modifying the planet's phase variation by the phase-dependent projected area of the ellipsoid, even for cases where the deformation is barely or not detectable. Atmospheric studies of UHJs through transmission spectroscopy could also be affected by tidal deformation. Lendl et al. (2017) suggested that planet deformation could be the cause of the enhanced features measured in the transmission spectra of WASP-103 b, but this remains to be understood and can be investigated using our ellipsoidal planet model.

**Acknowledgements.** This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. M.L. and B.A. acknowledge the financial support of the Swiss National Science Foundation under grant number PCEFP2_194576. S.C.C.B. acknowledges support from FCT -Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UIDB/04434/2020; UIDP/04434/2020, 2022.06962.PTDC.

## References

Akinsanmi, B., Barros, S. C., Santos, N. C., et al. 2019, A&A, 621, A117

Allard, F., Homeier, D., Freytag, B., et al. 2012, Phil. Trans. R. Soc. A: Math., Phys. Eng. Sci., 370, 2765

Arkhypov, O. V., Khodachenko, M. L., & Hanslmeier, A. 2021, A&A, 646, A136

Barros, S. C. C., Pollacco, D. L., Gibson, N. P., et al. 2012, MNRAS, 419, 1248

Barros, S. C. C., Akinsanmi, B., Boué, G., et al. 2022, A&A, 657, A52

Batygin, K., Bodenheimer, P., & Laughlin, G. 2009, ApJ, 704, L49

Baumeister, P., & Tosi, N. 2023, A&A, 676, A106

Baumeister, P., Padovan, S., Tosi, N., et al. 2020, ApJ, 889, 42

Bell, T. J., Zhang, M., Cubillos, P. E., et al. 2019, MNRAS, 489, 1995

Berardo, D., & de Wit, J. 2022a, ApJ, 935, 178

Berardo, D., & de Wit, J. 2022b, ApJ, 941, 155

Bourrier, V., Kitzmann, D., Kuntzer, T., et al. 2020, A&A, 637

Buhler, P. B., Knutson, H. A., Batygin, K., et al. 2016, ApJ, 821, 26

Burton, J. R., Watson, C. A., Fitzsimmons, A., et al. 2014, ApJ, 789, 113

Claret, A., & Southworth, J. 2022, A&A, 664, A128

Collins, K. A., Kielkopf, J. F., & Stassun, K. G. 2017, AJ, 153

Correia, A. C. M. 2014, A&A, 570, A5

Cowan, N. B., Machalek, P., Croll, B., et al. 2012, ApJ, 747, 82

Csizmadia, S., Smith, A. M. S., Kálmán, S., et al. 2023, A&A, 675, A106

Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e2019GL086572

Espinoza, N., & Jones, K. 2021, AJ, 162, 165

Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2013, ApJ, 772, 51

Hardy, R. A., Harrington, J., Hardin, M. R., et al. 2017, ApJ, 836, 143

Helled, R., Bodenheimer, P., Podolak, M., et al. 2014, in Protostars and Planets VI (University of Arizona Press)

Hellard, H., Csizmadia, S., Padovan, S., et al. 2019, ApJ, 878, 119

Hellard, H., Csizmadia, S., Padovan, S., Sohl, F., & Rauer, H. 2020, ApJ, 889, 66

Helling, C., Kawashima, Y., Graham, V., et al. 2020, A&A, 641, A178

Hestroffer, D. 1997, A&A, 327, 199

Hooton, M. J., de Mooij, E. J. W., Watson, C. A., et al. 2019, MNRAS, 486, 2397

Jones, K., & Espinoza, N. 2020, J. Open Source Softw., 7, 2382

Kellermann, C., Becker, A., & Redmer, R. 2018, A&A, 615, 39

Kramm, U., Nettelmann, N., Redmer, R., & Stevenson, D. J. 2011, A&A, 528, A18

Kramm, U., Nettelmann, N., Fortney, J. J., Neuhäuser, R., & Redmer, R. 2012, A&A, 538, A146

Kreidberg, L. 2018, in Handbook of Exoplanets (Cham: Springer International Publishing), 2083

Kreidberg, L., Line, M. R., Bean, J. L., et al. 2015, ApJ, 814, 66

Kreidberg, L., Line, M. R., Parmentier, V., et al. 2018, AJ, 156, 17

Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286

Leconte, J., Lai, D., & Chabrier, G. 2011, A&A, 528, A41

Lendl, M., Cubillos, P. E., Hagelberg, J., et al. 2017, A&A, 606

Loeb, A., & Gaudi, B. S. 2003, ApJ, 588, L117

Louden, T., & Kreidberg, L. 2018, MNRAS, 477, 2613

Love, A. E. H. 1911, Nature, 89, 471

Maxted, P. F. 2018, A&A, 616, 39

Mislis, D., Heller, R., Schmitt, J. H. M. M., & Hodgkin, S. 2012, A&A, 538, A4

Morello, G., Tsiaras, A., Howarth, I. D., & Homeier, D. 2017, AJ, 154, 111

Morris, S. L. 1985, ApJ, 295, 143

Newville, M., Otten, R., Nelson, A., et al. 2020, https://zenodo.org/records/3814709

Padovan, S., Spohn, T., Baumeister, P., et al. 2018, A&A, 620, A178

Parviainen, H. 2018, in Handbook of Exoplanets (Cham: Springer International Publishing), 1567

Parviainen, H., Wilson, T. G., Lendl, M., et al. 2022, A&A, 668, A93

Pluriel, W., Zingales, T., Leconte, J., & Parmentier, V. 2020, A&A, 636, A66

Ragozzine, D., & Wolf, A. S. 2009, ApJ, 698, 1778

Rauer, H., & Heras, A. M. 2018, in Handbook of Exoplanets (Cham: Springer International Publishing), 1309

Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249

Russell, H. N. 1912, ApJ, 35, 315

Sabadini, R., & Vermeersen, B. 2004, Global Dynamics of the Earth: Applications of Normal Mode Relaxation Theory to Solid-Earth Geophysics (Kluwer Academic Publishers)

Shporer, A. 2017, PASP, 129, 072001

Shporer, A., Wong, I., Huang, C. X., et al. 2019, AJ, 157, 178

Speagle, J. S. 2020, MNRAS, 493, 3132

Stevenson, K. B., Bean, J. L., Seifahrt, A., et al. 2014, AJ, 147, 161

Wahl, S. M., Thorngren, D., Lu, T., & Militzer, B. 2021, ApJ, 921, 105

Wong, I., Kitzmann, D., Shporer, A., et al. 2021, AJ, 162, 127

Wong, I., Shporer, A., Vissapragada, S., et al. 2022, AJ, 163, 175

Zhang, X., & Showman, A. P. 2017, ApJ, 836, 73

Zieba, S., Zilinskas, M., Kreidberg, L., et al. 2022, A&A, 664, A79
**Appendix A: Additional figures**

<figure><img src="image_10.png" /><figcaption>**Fig. A.1.** Comparison of the numeric and analytic computation of the planetary phase variation signal. The spherical planet signal is equivalent to a cosine function.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. A.2.** Same as Fig. 5, but for HST/WFC3/ simulations.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. A.3.** Same as Fig. 6, but for simulated 200 ppm/2 min noise observations for TESS.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. A.4.** Same as Fig. 6, but for simulated 30 ppm/2 min noise observations of PLATO.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. A.5.** Same as Fig. 6, but for simulated 85 ppm/2 min noise observations for CHEOPS.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. A.6.** Same as Fig. 6, but for simulated 110 ppm/2 min noise observations for HST/WFC3-G141.</figcaption></figure>