Signatures of Obliquity in Thermal Phase Curves of Hot Jupiters

TL;DR

This paper investigates how large planetary obliquities influence thermal phase curves of hot Jupiters, proposing indirect detection methods and analyzing specific exoplanets for signs of high obliquity.

Contribution

It introduces a thermal radiative model to identify obliquity signatures in phase curves and applies it to known hot Jupiters to constrain their obliquities.

Findings

CoRoT-2 b shows potential high obliquity signatures.

No significant obliquity evidence for HD 149026 b.

WASP-12 b's phase curve is too complex for clear conclusions.

Abstract

Recent work suggests that many short-period extrasolar planets may have spin obliquities that are significantly tilted with respect to their orbital planes. These large obliquities are a natural outcome of ``secular spin-orbit resonance'', a configuration in which the planetary spin precession frequency matches the frequency of orbit nodal regression, or a Fourier component thereof. While exoplanet spin obliquities have not yet been measured directly, they may be detectable indirectly through their signatures in various observations, such as photometric measurements across the full phase of a planet's orbit. In this work, we employ a thermal radiative model to explore how large polar tilts affect full-phase light curves, and we discuss the range of unique signatures that are expected to result. We show that the well-studied short-period planets HD 149026 b, WASP-12 b, and CoRoT-2 b all exhibit phase curve features that may arise from being in high-obliquity states. We also constrain the parameters and assess the detectability of hypothetical perturbing planets that could maintain the planets in these states. Among the three planets considered, CoRoT-2 b has the tightest constraints on its proposed obliquity ($45.8^{\circ} \pm 1.4^{\circ}$) and axial orientation. For HD 149026 b, we find no significant evidence for a non-zero obliquity, and the phase curve of WASP-12 b is too complicated by strong tidal distortions for a conclusive assessment.

Tables

Table 1. Table 1: Planetary System Properties

	HD 149026 b	WASP-12 b	CoRoT-2 b
$P$ (days)	$2.8758911\pm 2.5\times {10}^{-6}$	$1.09142119\pm 2.1\times {10}^{-7}$	$1.7429935\pm 1.0\times {10}^{-6}$
$e$	$0$	$0.0447\pm 0.0043$	${0.0143}_{-0.0076}^{+0.0077}$
$\varpi ({}_{}{}^{\circ })$	N/A	${272.7}_{-1.3}^{+2.4}$	${102}_{-5}^{+17}$
$M_P$ $\left(M_J\right)$	${0.368}_{-0.014}^{+0.013}$	$1.43\pm 0.14$	$3.47\pm 0.22$
$R_P$ $\left(R_J\right)$	${0.813}_{-0.025}^{+0.027}$	$1.825\pm 0.094$	${1.466}_{-0.044}^{+0.042}$
$M_{\star }$ $\left(M_{\odot }\right)$	$1.345\pm 0.020$	$1.280\pm 0.05$	$0.96\pm 0.08$
$R_{\star }$ $\left(R_{\odot }\right)$	${1.541}_{-0.042}^{+0.046}$	$1.630\pm 0.08$	${0.906}_{-0.027}^{+0.026}$
$T_{\text{eff}}$ (K)	$6160\pm 50$	${6300}_{-100}^{+200}$	$5625\pm 120$
Ref.	(1)	(2)–(4)	(5)

Table 2. Table 2: Best-Fit Parameters from Radiative Model

	HD 149026 b		WASP-12 b		CoRoT-2 b
Parameter	3.6 $\mu$m	4.5 $\mu$m	3.6 $\mu$m	4.5 $\mu$m	4.5 $\mu$m
Sub-synchronous Rotation
${\tau }_{\mathrm{rad}}$ (hr)	${7.7}_{-4.6}^{+10.8}$	${2.4}_{-1.3}^{+24.0}$	${0.01}_{-0.01}^{+0.01}$	${65.5}_{-0.3}^{+24.3}$	${23.9}_{-3.1}^{+2.5}$
$T_0$ (K)	${877}_{-631}^{+388}$	${1388}_{-158}^{+155}$	${427}_{-N/A}^{+166}$	${2458}_{-141}^{+84}$	${92}_{-52}^{+209}$
$A$		${0.64}_{-0.05}^{+0.09}$		${0.36}_{-0.35}^{+0.09}$
$ϵ$ (∘)	${93.9}_{-32.5}^{+45.1}$	${4.2}_{-1.9}^{+47.7}$	${87.9}_{-8.5}^{+0.6}$	${91.2}_{-3.1}^{+70.1}$	$45.8\pm 1.4$
$\zeta$ (∘)	${10.6}_{-170.42}^{+143.7}$	$-{29.0}_{-149.2}^{+113.0}$	$-{39.4}_{-0.9}^{+4.4}$	$-{73.0}_{-58.6}^{+21.2}$	$-{82.5}_{-7.1}^{+10.8}$
Free Rotation
$P_{\mathrm{rot}}/P_{\mathrm{orb}}$	$1.09\pm 0.06$	${0.39}_{-0.21}^{+0.08}$	${0.77}_{-0.04}^{+0.14}$	$0.94\pm 0.01$	${1.13}_{-0.02}^{+0.04}$
${\tau }_{\mathrm{rad}}$ (hr)	${97.2}_{-0.2}^{+5.1}$	${3.7}_{-2.2}^{+0.4}$	${12.7}_{-2.7}^{+32.2}$	${40.0}_{-3.4}^{+8.1}$	${16.9}_{-4.4}^{+4.7}$
$T_0$ (K)	${1172}_{-111}^{+334}$	${1223}_{-894}^{+165}$	${1567}_{-63}^{+134}$	${2204}_{-14}^{+40}$	${102}_{-75}^{+219}$
$A$	${0.01}_{-0.01}^{+0.39}$	${0.44}_{-0.34}^{+0.05}$	${0.09}_{-0.08}^{+0.02}$
$ϵ$ (∘)	${4.2}_{-2.7}^{+57.2}$	${36.4}_{-3.4}^{+28.3}$	${0.8}_{-0.4}^{+6.0}$	${21.1}_{-10.1}^{+8.5}$	${2.0}_{-1.2}^{+8.1}$
$\zeta$ (∘)	$-{4.7}_{-30.2}^{+7.3}$	${69.8}_{-29.5}^{+3.0}$	$-{67.0}_{-N/A}^{+123.4}$	${46.6}_{-12.1}^{+58.3}$	$-{32.3}_{-20.8}^{+47.1}$

Table 3. Table 3: Model Comparison by Number of Free Parameters

Obliquity		Rotation		$\Delta$AIC
Case	$ϵ$ (∘)	Case	$P_{\mathrm{rot}}/P_{\mathrm{PSR}}$	3.6 $\mu$m	4.5 $\mu$m	Combined
HD 149026 b
Fixed zero	$0^{\circ }$	Synchronous	1	2.62	0	1.37
Fixed zero	$0^{\circ }$	Free	$1.11\pm 0.05$, ${0.51}_{-0.03}^{+0.28}$	0	1.26	0
$\lambda$-consistent	${89.3}_{-31.8}^{+30.2}$	Sub-synchronous	43.41 ($>1.20$)	1.80	4.11	0.51
$\lambda$-consistent	${119.0}_{-49.1}^{+41.5}$	Free	${0.37}_{-0.04}^{+2.25}$, ${1.21}_{-0.68}^{+0.40}$	3.51	6.32	4.39
Free	${93.9}_{-32.5}^{+45.1}$, ${4.2}_{-1.9}^{+47.7}$	Sub-synchronous	$7.31$ ($>1.04$), ${1.00}_{-0.00}^{+0.12}$	1.78	4.12	4.63
Free	${4.2}_{-2.7}^{+57.2}$, ${36.4}_{-3.4}^{+28.3}$	Free	$1.09\pm 0.06$, ${0.39}_{-0.21}^{+0.08}$	4.09	5.35	8.13
WASP-12 b
Fixed zero	$0^{\circ }$	Synchronous	1	593.24	91.18	680.13
Fixed zero	$0^{\circ }$	Free	${0.91}_{-0.00}^{+0.01}$, $0.95\pm 0.01$	0	6.16	1.87
$\lambda$-consistent	${19.4}_{-1.6}^{+1.7}$	Sub-synchronous	$1.002\pm 0.001$	592.90	81.67	666.17
$\lambda$-consistent	${11.1}_{-7.9}^{+4.5}$	Free	${0.910}_{-0.001}^{+0.009}$, ${0.945}_{-0.005}^{+0.010}$	68.37	6.63	66.57
Free	${93.9}_{-32.5}^{+45.1}$, ${4.2}_{-1.9}^{+47.7}$	Sub-synchronous	${1.0015}_{-0.0010}^{+0.0005}$, ${1.0039}_{-0.0009}^{+0.0011}$	586.68	76.60	658.98
Free	${4.4}_{-2.2}^{+3.3}$, ${27.4}_{-10.8}^{+6.2}$	Free	$0.91\pm 0.004$, $0.95\pm 0.01$	4.31	0	0
CoRoT-2 b
Fixed zero	$0^{\circ }$	Synchronous	1		259.41
Fixed zero	$0^{\circ }$	Free	${1.10}_{-0.03}^{+0.02}$		0
Free	$45.8\pm 1.4$	Sub-synchronous	$1.07\pm 0.01$		7.56
Free	${2.0}_{-1.2}^{+8.1}$	Free	${1.13}_{-0.02}^{+0.04}$		4.06

Table 4. Table 4: Orbital and Transit Properties for Kepler-107

Planet	$P$ (days)	$a$ (AU)	$e$	$M_p$ ($M_J$)	$R_p$ ($R_J$)	$b$	$T_{\mathrm{tran}}$ ($\mathrm{BJD}-2454900$)
b	$3.179997\pm 1.1\times {10}^{-5}$	0.044	${0.020}_{-0.020}^{+0.200}$	0.01167	$0.139\pm 0.005$	$0.34\pm 0.24$	$66.49904\pm 0.00200$
c	$4.901425\pm 1.6\times {10}^{-5}$	0.059	${0.020}_{-0.020}^{+0.260}$	0.0133	$0.161\pm 0.016$	$0.78\pm 0.29$	$71.60742\pm 0.00174$
d	$7.958203\pm 1.04\times {10}^{-4}$	0.082	${0.14}_{-0.14}^{+0.25}$	0.00371	$0.095\pm 0.005$	$0.27\pm 0.24$	$70.79968\pm 0.00612$
e	$14.749049\pm 3.4\times {10}^{-5}$	0.123	${0.020}_{-0.020}^{+0.180}$	0.0360	$0.308\pm 0.022$	$0.90\pm 0.28$	$71.77998\pm 0.00134$

Table 5. Table 5: System parameters used in developing the Figure 6 constraints on the perturbing planet parameters. The quantities k 2 ⁣ ⋆ subscript 𝑘 2 ⋆ k_{2\star} , k 2 subscript 𝑘 2 k_{2} , and C 𝐶 C are estimates. α syn subscript 𝛼 syn \alpha_{\mathrm{syn}} is the spin axis precession constant that would result from the aggregate of system parameters in the case of synchronous rotation ( ω = α syn 𝜔 subscript 𝛼 syn \omega=\alpha_{\mathrm{syn}} ). g ⋆ subscript 𝑔 ⋆ g_{\star} is the approximate frequency of nodal recession induced by stellar oblateness. Note that α syn ≫ g ⋆ much-greater-than subscript 𝛼 syn subscript 𝑔 ⋆ \alpha_{\mathrm{syn}}\gg g_{\star} in all cases.

	HD 149026 b	WASP-12 b	CoRoT-2 b
$M_{\star }\left[M_{\odot }\right]$	${1.345}^{\left(1\right)}$	${1.36}^{\left(2\right)}$	${0.97}^{\left(3\right)}$
$R_{\star }\left[R_{\odot }\right]$	${1.541}^{\left(1\right)}$	${1.63}^{\left(2\right)}$	${0.90}^{\left(3\right)}$
$P_{\star }\left[\mathrm{days}\right]$	$\sim {13}^{\left(4^{\ast }\right)}$	${36}^{\left(5\right)}$	${4.5}^{\left(3\right)}$
$k_{2\star }$	0.01(6)	0.01(6)	0.01(6)
$a\left[\mathrm{AU}\right]$	${0.042}^{\left(4\right)}$	${0.02299}^{\left(7\right)}$	${0.0281}^{\left(3\right)}$
$M_p\left[M_{\mathrm{Jup}}\right]$	${0.36}^{\left(4\right)}$	${1.41}^{\left(7\right)}$	${3.3}^{\left(3\right)}$
$R_p\left[R_{\mathrm{Jup}}\right]$	${0.725}^{\left(4\right)}$	${1.89}^{\left(2\right)}$	${1.466}^{\left(3\right)}$
$k_2$	0.1	0.1	0.1
$C$	0.2	0.2	0.2
$2\pi /{\alpha }_{\mathrm{syn}}\left[\mathrm{yr}\right]$	14.3	0.2	4.0
$2\pi /g_{\star }\left[\mathrm{yr}\right]$	$2.4\times {10}^5$	$1.7\times {10}^5$	$8.7\times {10}^4$

Table 6. Table 6: RV semi-amplitudes of the known hot Jupiters in each system, along with RV jitters of the fits.

	$K_1$ [m/s]	${\sigma }_{\mathrm{jit}}$ [m/s]	Reference
HD 149026 b	$38.5\pm 1.2$	$5.8\pm 0.6$	(1)
WASP-12 b	${219.9}_{-2.1}^{+2.2}$	${9.1}_{-1.3}^{+1.8}$	(2)
CoRoT-2 b	${568}_{-22}^{+23}$	${40}_{-10}^{+14}$	(2)

Equations

$$r\left(t\right)=a\left(\frac{1-e^2}{1+e\cos \nu }\right)$$

$$M\left(t\right)=E-e\sin E.$$

$$\begin{array}{rl}\cos \nu & =\frac{\cos E-e}{1-e\cos E}\\ \sin \nu & =\frac{\sqrt{1-e^2}\sin E}{1-e\cos E}.\end{array}$$

$$T_{\text{eq}}^4\left(\varphi ,\theta ,t\right)=\left(1-A\right)\left(\frac{L_{\star }}{4\pi \sigma r^2}\right)\cos {\alpha }_{\star }+T_0^4$$

$$\cos {\alpha }_{\star }=\{\begin{array}{ll}\hat{n}\cdot {\hat{r}}_{\star },& \hat{n}\cdot {\hat{r}}_{\star }\ge 0\\ 0,& \hat{n}\cdot {\hat{r}}_{\star }<0\end{array}$$

$$\displaystyle\begin{split}\phi_{\star}=\ &\phi_{\star}\!\left(t_{0}\right)-\Big{\{}\omega_{\mathrm{rot}}\left(t-t_{0}\right)+\\ &\arctan\left\{\cos\epsilon\tan\left[\left(\nu-\zeta\right)-\nu\!\left(t_{0}\right)\right]\right\}\Big{\}}\\ \theta_{\star}=\ &\sin^{-1}\left[\sin\epsilon\cos\left(\nu-\zeta\right)\right].\end{split}$$

$$\dot{T}\left(\varphi ,\theta ,t\right)=\frac{T_{\text{eq}}}{4{\tau }_{\text{rad}}}\left\{1-{\left[\frac{T\left(\varphi ,\theta ,t\right)}{T_{\text{eq}}}\right]}^4\right\},$$

$$\bar{F}\left(t\right)=\frac{1}{\pi }{\left(\frac{R_{\mathrm{p}}}{R_{\star }}\right)}^2\frac{\iiint w{B}_{\lambda }\left(T\right)Vd\lambda d\theta d\varphi }{\int w{B}_{\lambda }\left(T_{\star }\right)d\lambda }$$

$$V=\{\begin{array}{ll}\hat{n}\cdot {\hat{r}}_{\mathrm{obs}},& \hat{n}\cdot {\hat{r}}_{\mathrm{obs}}\ge 0\\ 0,& \hat{n}\cdot {\hat{r}}_{\mathrm{obs}}<0\end{array}$$

$$\begin{array}{rl}{\varphi }_{\mathrm{obs}}& ={\varphi }_{\star }\left(t_{\mathrm{ecl}}\right)-{\omega }_{\mathrm{rot}}\left(t-t_{\mathrm{ecl}}\right)\\ {\theta }_{\mathrm{obs}}& ={\sin }^{-1}\left[\sin ϵ\cos \left({\nu }_{\mathrm{ecl}}-\zeta \right)\right]\end{array}$$

$$\displaystyle\begin{split}\Phi\!\left(\vec{r}\right)=-\frac{GM_{\star}}{a}\Bigg{\{}&\left[\left(\frac{z}{a}+\frac{\xi}{1+\xi}\right)^{2}+\left(\frac{x^{2}+y^{2}}{a^{2}}\right)\right]^{-1/2}\\ +&\left[\left(\frac{z}{a}+\frac{1}{1+\xi}\right)^{2}+\left(\frac{x^{2}+y^{2}}{a^{2}}\right)\right]^{-1/2}\\ +&\frac{1+\xi}{2}\left(\frac{x^{2}+z^{2}}{a^{2}}\right)\Bigg{\}}\end{split}$$

$$\displaystyle\begin{split}S_{ij}&=\,B^{2}_{\mathrm{p}}\int_{\phi_{i}}^{\phi_{j}}\!\int_{\theta_{i}}^{\theta_{j}}\!\cos\theta\,\Big{\{}1\\ &+\left(\chi^{2}-1\right)\left[\frac{1}{4}f\!\left(\phi,\theta\right)+2\left(\hat{r}\cdot\hat{r}_{\ell}\right)\right]\\ &+\left(\chi^{2}-1\right)^{2}\left(\hat{r}\cdot\hat{r}_{\ell}\right)^{2}\Big{\}}^{1/2}\,d\theta\,d\phi,\end{split}$$

$$\begin{array}{rl}f\left(\varphi ,\theta \right)& ={\left\{\sin \theta \cos {\theta }_{\ell }\left[\cos \left(\varphi -{\varphi }_{\ell }\right)-1\right]-\sin \left(\theta -{\theta }_{\ell }\right)\right\}}^2\\ & +{\cos }^2{\theta }_{\ell }\left[\sin \left(\varphi -{\varphi }_{\ell }\right)+1\right].\end{array}$$

$$f\left(x,y,z\right)=\frac{x^2+y^2}{B_{\mathrm{p}}^2}+\frac{z^2}{A_{\mathrm{p}}^2}.$$

$$\begin{array}{rl}{\varphi }_{\ell }& ={\varphi }_{\star }\left(\nu \to \nu -\lambda \right)\\ {\theta }_{\ell }& ={\theta }_{\star }\left(\nu \to \nu -\lambda \right).\end{array}$$

$$\begin{array}{rl}\hat{x}& =\frac{{\hat{r}}_{\star }-\cos \lambda {\hat{r}}_{\ell }}{\sin \lambda }\\ \hat{y}& =\frac{{\hat{r}}_{\ell }\times {\hat{r}}_{\star }}{\sin \lambda }\\ \hat{z}& ={\hat{r}}_{\ell }.\end{array}$$

$$\begin{array}{rl}\hat{n}& =\vec{\nabla }f/\Vert \vec{\nabla }f\Vert \\ & ={\left\{1+g^{-1}\left[\chi ,\left(\hat{r}\cdot {\hat{r}}_{\ell }\right)\right]\right\}}^{-1/2}\left[\left(\hat{r}\cdot \hat{x}\right)\hat{x}+\left(\hat{r}\cdot \hat{y}\right)\hat{y}\right]\\ & +{\left\{1+g\left[\chi ,\left(\hat{r}\cdot {\hat{r}}_{\ell }\right)\right]\right\}}^{-1/2}\hat{z}\end{array}$$

$$g\left[\chi ,\left(\hat{r}\cdot {\hat{r}}_{\ell }\right)\right]\equiv {\chi }^2\left[\frac{1-{\left(\hat{r}\cdot {\hat{r}}_{\ell }\right)}^2}{{\left(\hat{r}\cdot {\hat{r}}_{\ell }\right)}^2}\right]$$

$$\begin{array}{rl}\cos {\alpha }_{\ell }& =\hat{n}\cdot {\hat{r}}_{\star }\\ & =\sin \lambda {\left(1+g^{-1}\right)}^{-1/2}\left[\hat{r}\cdot {\hat{r}}_{\star }-\cos \lambda \left(\hat{r}\cdot {\hat{r}}_{\ell }\right)\right]\\ & +\cos \lambda {\left(1+g\right)}^{-1/2}\left(\hat{r}\cdot {\hat{r}}_{\ell }\right)\end{array}$$

$$\begin{array}{rl}V& =\hat{n}\cdot {\hat{r}}_{\mathrm{obs}}\\ & =\frac{\cos \left(\nu -{\nu }_{\mathrm{ecl}}\right)-\cos \left[\left(\nu -\lambda \right)-{\nu }_{\mathrm{ecl}}\right]}{\sin \lambda }{\left(1+g^{-1}\right)}^{-1/2}\\ & \times \left[\hat{r}\cdot {\hat{r}}_{\star }-\cos \lambda \left(\hat{r}\cdot {\hat{r}}_{\ell }\right)\right]\\ & +\cos \left[\left(\nu -\lambda \right)-{\nu }_{\mathrm{ecl}}\right]{\left(1+g\right)}^{-1/2}\left(\hat{r}\cdot {\hat{r}}_{\ell }\right).\end{array}$$

$${\omega }_{\mathrm{eq}}=n\frac{2\cos ϵ}{1+{\cos }^2ϵ}$$

$$\mathrm{AICc}\equiv 2\left[k-\ln \mathcal{L}+\frac{k\left(k+1\right)}{n-k-1}\right],$$

$$|\sin \Omega |\ge b\left[\frac{R_{\star }+R_{\mathrm{P}}}{r_{\mathrm{tr}}\left(e,\varpi \right)}\right]$$

$$r_{\mathrm{tr}}\equiv a\left(\frac{1-e^2}{1+e\sin \varpi }\right)$$

$$\sin i|\sin \Omega |=b\left[\frac{R_{\star }+R_{\mathrm{P}}}{r_{\mathrm{tr}}\left(e,\varpi \right)}\right].$$

$$\alpha =\frac{1}{2}\frac{M_{\star }}{M_{\mathrm{p}}}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^3\frac{k_2}{C}\omega .$$

$$J_2=\frac{{\omega }^2{R_{\mathrm{p}}}^3}{3G{M}_{\mathrm{p}}}{k}_2.$$

$$\begin{array}{rl}{g}_{{}_{\mathrm{LL}}}=& -\frac{1}{4}{b}_{3/2}^{(1)}({\alpha }_{12}){\alpha }_{12}\times \\ & \left(n_1\frac{M_{p2}}{M_{\star }+M_{p1}}{\alpha }_{12}+n_2\frac{M_{p1}}{M_{\star }+M_{p2}}\right),\end{array}$$

$$b_{3/2}^{(1)}({\alpha }_{12})=\frac{1}{\pi }{\int }_0^{2\pi }\frac{\cos \psi }{(1-2\alpha \cos \psi +{\alpha }^2{)}^{3/2}}d\psi .$$

$$g_{\star }=n\frac{k_{2\star }}{2}{\left(\frac{{\omega }_{\star }}{n}\right)}^2{\left(\frac{R_{\star }}{a}\right)}^5.$$

Full text

methods: data analysis, methods: numerical, planets and satellites: atmospheres, planets and satellites: dynamical evolution and stability, planets and satellites: individual (HD 149026 b, WASP-12 b, CoRoT-2 b), techniques: photometric

††journal: AJ††software: Astropy (Astropy Collaboration et al. 2013), Colorpy (Kness 2008), python-colormath (Taylor 2018), Jupyter (Kluyver et al. 2016), Matplotlib (Hunter 2007), Numpy (van der Walt et al. 2011), Paletton (Staníček 2018), REBOUND (Rein & Liu 2012; Rein & Spiegel 2015), Scipy (Jones et al. 2001), WebPlotDigitizer (Rohatgi 2017)

1 Introduction

Analysis of the full-phase photometry of transiting extrasolar planets has generated a number of insights. By tracking the emission and reflection of light from planets as they trace through their orbits, one may derive important clues regarding the planets’ atmospheric compositions, surface flow patterns, atmospheric thermal responses, cloud coverage, day-to-night heat redistribution efficiency, and more. High signal-to-noise detections of full-phase photometry are now routinely performed in both the optical and in the near-infrared, and the technique has been employed for both detection and characterization (see, e.g. Heng & Showman 2015; Shporer 2017; Deming & Seager 2017).

The climate properties probed by full-phase light curves are highly sensitive to the planet’s orbital and spin geometry. For instance, many short-period ($P\lesssim 5\,{\rm d}$) planets are subject to strong tidal evolution that has produced complete or near-complete orbital circularization (Rodríguez & Ferraz-Mello 2010). Moreover, because time scales for planetary tidal spin evolution are generally shorter than for orbital evolution, it is usually assumed that short-period planets have zero obliquity, and are spinning synchronously for orbits with $e=0$ (see e.g. Gladman et al. 1996) or pseudo-synchronously for orbits with finite eccentricity (see e.g. Hut 1981; Ivanov & Papaloizou 2007).

Time-resolved photometry of assumed-synchronous planets on circular orbits suggests that the peak infrared emission from the planet usually lies eastward of the sub-stellar point (e.g. Knutson et al. 2007, 2009a; Cowan & Agol 2011; Knutson et al. 2012; Cowan et al. 2012b; Zellem et al. 2014; Schwartz et al. 2017; Zhang et al. 2018). By contrast, phase curves in the optical tend to suggest that peak reflectivity occurs westward of the sub-stellar point and closer to the morning terminator (Shporer 2017). These observations are generally interpreted to imply eastward circumplanetary flow and cloud decks that burn off when advected into the direct beam of instellation.

This composite picture can be tested when the phase curves of eccentric orbits are tracked (e.g. Laughlin et al. 2009; Adams et al. 2019). A planet with a non-circular orbit cannot be fully tidally de-spun; the surface flows are thus dynamically responsive. Bulk properties of the atmosphere such as the radiative response timescale can be directly inferred (e.g. Cowan et al. 2012b; de Wit et al. 2016; Adams & Laughlin 2018).

Just as the eccentricity produces observable signatures in the full-phase photometry, so too will the presence of a significant obliquity. It is thus of interest to examine possible obliquity-induced signatures in the context of the data sets that are currently available.

Several authors have developed a mathematical formalism to predict optical observational effects of the relationship between a planet’s spin, orbit, and viewing geometry. Kawahara & Fujii (2010) constructed a geometric framework for mapping planets (including oblique ones) in reflected/scattered light. This framework was later extended to account for effects such as cloud cover (Kawahara & Fujii 2011) and generally inclined orbits (Fujii & Kawahara 2012). Signatures of obliquity have been predicted in Fourier analyses of photometry for directly-imaged planets (Kawahara 2016). Schwartz et al. (2016) demonstrated the ability to infer spin axis orientation for general albedo maps, and outlined a feasible minimum observing baseline for making robust inferences of the spin axis orientation. Farr et al. (2018) introduced exocartographer111https://github.com/bfarr/exocartographer, a software package that generates reflected light photometry for an arbitrary albedo map and spin geometry and fits a variety of time-sampled data. Another software package, STARRY222https://rodluger.github.io/starry/ (Luger et al. 2019), provides computationally efficient determinations of full-phase light curves, occultation signatures, and transit signals by using a global planetary surface pattern expressed as a sum of spherical harmonics.

Here we consider a similar geometric framework in the near-infrared, where the thermal emission from the planet — rather than reflected light — should be the primary source of the observed flux. Cowan et al. (2012b) considered the effects of eccentricity and obliquity for the phase photometry of an Earth-like planet and demonstrated that one can infer the characteristic thermal time scales. They also concluded that a combination of optical and infrared observations would be necessary to accurately measure bulk atmospheric conditions. Most recently, Ohno & Zhang (2019a, b) developed a comprehensive shallow water model for describing the atmospheric dynamics and resulting thermal phase variations of planets with arbitrary spin period, spin obliquity, and orbital eccentricity. In short, because these three properties all determine the instantaneous sub-stellar location, the resulting thermal phase variations will be shaped by their combined influence.

While one can model thermal variations for any choice of orbit and spin orientation, we will first restrict ourselves to scenarios that are physically feasible for systems with close-in giant planets. In addition to representing a dynamically plausible configuration, accurate phase curves for short-period giant planets already exist. Moreover, if we infer a particular spin-orbit geometry for a planet from its observed photometry, we may place constraints on the system architecture and make potentially observable dynamical predictions.

In this work, we focus on three short-period giant planets: HD 149026 b, WASP-12 b, and CoRoT-2 b. These all have full-phase thermal light curves obtained with Spitzer, and they have features (discussed in detail for HD 149026 b and WASP-12 b in Adams & Laughlin 2018) that suggest they are potential high obliquity candidates. HD 149026 b, WASP-12 b, and CoRoT-2 b do not currently have any known planetary companions. They would therefore require special dynamical states to have anything other than circular orbits, spin-orbit synchronization, and zero obliquities.

One such dynamical state that may maintain a large obliquity is a “Cassini state”, an equilibrium configuration where the planet’s spin vector stays fixed in the reference frame of its precessing orbit (Colombo 1966; Peale 1969; Ward 1975). As an instance of a Cassini state, a secular spin-orbit resonance involves an average commensurability between the frequency of the planet’s spin axis precession and the frequency of its orbit nodal recession (or a Fourier component thereof). The orbital recession may be provided by a number of sources, such as another planet in the system, the stellar quadrupolar gravitational potential (Fabrycky et al. 2007; Millholland & Laughlin 2019), or, early on in the system’s lifetime, the protoplanetary disk (Millholland & Batygin 2019). The result is a stable state in which a planet may maintain a non-synchronous spin and non-zero obliquity, even in the presence of strong tides.

Cassini states have been invoked to explain, for example, the co-precession of the lunar spin vector and the lunar orbit normal (Delaunay 1860), the obliquities of Saturn (Ward & Hamilton 2004; Hamilton & Ward 2004) and Jupiter (Ward & Canup 2006), and the spin precession state of Mercury (Peale 2006). They have also been considered as a mechanism to inflate the radii of hot Jupiters (Winn & Holman 2005), albeit with disputed feasibility (Levrard et al. 2007; Fabrycky et al. 2007). Recently, Millholland & Laughlin (2019) showed that high-obliquity Cassini states might be common for planets in short-period, compact systems. In this paper, we examine whether HD 149026 b, WASP-12 b, and CoRoT-2 b may have large obliquities due to their participation in secular spin-orbit resonances, and we examine the potential signatures of such states in their thermal full-phase light curves.

This paper is organized as follows. In §2, we introduce three unusual planets – HD 149026 b, WASP-12 b, and CoRoT-2 b – and discuss why they are viable candidates for obliquity investigations. In §3 we review these planets’ near-infrared photometric observations, including the analyses of their light curve morphologies. §4 describes the thermal model, which accounts for a general rotation rate and spin axis orientation. This framework is then employed in §5 to re-analyze the Spitzer phase photometry. In order to build a physical framework where we could plausibly observe a single-transiting, close-in, oblique planet, we make a feasibility assessment in §6 wherein a compact, nearly co-planar multi-planet system evolves on short timescales to states of high mutual inclination, with only one planet capable of transiting. Finally, §7 provides a dynamical analysis of the three case study systems and explores the possibility that a Cassini state may exist between the known planets and an as-yet undiscovered, non-transiting companion.

2 Three Unusual Planets

HD 149026 b orbits its subgiant host in 2.88 days; the planet is slightly more massive (0.38 MJ) but smaller in size (0.74 RJ) than Saturn (Stassun et al. 2017). Its formation is still a matter of ongoing study, since the bulk density is quite high compared with other short-period giant planets. Sato et al. (2005), who first measured the planet’s radius during transit, proposed that a high core mass of $\sim\!67$ M⊕ could explain the measurements. Subsequent assessments have inferred similarly high core masses in the range of $\sim\!50$–110 M⊕ (Fortney et al. 2006; Ikoma et al. 2006; Broeg & Wuchterl 2007; Burrows et al. 2007) using combinations of atmospheric and interior modeling. Ikoma et al. (2006) proposed that HD 149026 b’s high metallicity and modest H/He envelope might either be explained via planetesimal capture and a limited gas supply, if its current state was obtained prior to disk dissipation, or a combination of envelope photoevaporation, Roche lobe overflow, and major collisions after disk dissipation. Zhang et al. (2018) find that the best-fit phase offsets from the Spitzer 3.6 $\mu$m and 4.5 $\mu$m full phase photometry are both significantly different from zero and in disagreement with each other. While Zhang et al. (2018) suggested uncharacterized instrumental systematics as a potential cause of the disparity, it is also worth investigating whether unusual global atmospheric dynamics may be at play.

WASP-12 b is another short-period ($1.09$ days) giant planet, which is tidally distorted due to its density and proximity to its host star (Li et al. 2010; Lai et al. 2010). There is spectroscopic evidence that the planet is overflowing its Roche lobe (Fossati et al. 2010; Haswell et al. 2012; Fossati et al. 2013; Jackson et al. 2017), and hydrodynamic simulations (Debrecht et al. 2018) suggest that it is undergoing significant atmospheric mass loss. von Essen et al. (2019) have recently analyzed time variability in the measured optical eclipse depths; with neither cloud albedo nor high temperatures tenable to explain the variability without invoking extreme physical magnitudes, they point cautiously to additional occultation from the atmospheric mass loss as a possible mechanism. Extensive transit observations have also been made (detailed further in §3), whose rapidly advancing ephemerides suggest either apsidal precession of an eccentric orbit or orbital decay (Maciejewski et al. 2016; Patra et al. 2017). Recent observations (Maciejewski et al. 2018) and theoretical investigations (Bailey & Goodman 2019) are in favor of the orbital in-spiral scenario. Millholland & Laughlin (2018) proposed that this rapid orbital decay may be due to tidal dissipation in the planet that is strongly enhanced by a high obliquity state.

CoRoT-2 b is a $3.3\ M_{\mathrm{J}}$ hot Jupiter with a 1.74-day orbital period and $1.47\ R_{\mathrm{J}}$ radius that has consistently been measured as anomalously inflated333It is important to note that the host star is quite active (Lanza et al. 2009; Huber et al. 2009), and this activity can bias measurements of the planetary radius from transit observations (Czesla et al. 2009; Huber et al. 2010; Bruno et al. 2016). (Alonso et al. 2008; Gillon et al. 2010; Southworth 2011). Previous proposed explanations include that the system is very young ($\sim\!30$–40 Myr) and the planet has not yet fully contracted gravitationally, or it is a bit older ($\sim\!130$–500 Myr) and the radius inflation is a transient response from a recent collision (Guillot & Havel 2011). Alternatively, the extreme radius inflation may be the result of obliquity tides, a possibility we discuss later in §8. In addition, Dang et al. (2018) reported Spitzer thermal phase observations that robustly indicate that the day-side hotspot is westward of the sub-stellar point, in contrast to the eastward offset that planets typically exhibit. Magnetohydrodynamic effects may be one mechanism for generating an unusual westward offset, but recently Hindle et al. (2019) inferred that such effects are very unlikely to be a viable mechanism for the specific westward offsets seen at 4.5 $\mu$m in CoRoT-2 b. Dang et al. (2018) considered that it may arise from non-synchronous rotation, which, as discussed in §7, would be a consequence of a high-obliquity state.

3 Photometry

HD 149026 b has been observed in transit in Strömgren $b$ and $y$ (Sato et al. 2005; Winn et al. 2008), $g$ and $r$ (Charbonneau et al. 2006), NICMOS (1.1–2.0 $\mu$m) on the Hubble Space Telescope (Carter et al. 2009), and the 8.0 $\mu$m channel of the Infrared Array Camera (IRAC) on the Spitzer Space Telescope (Nutzman et al. 2009). Secondary eclipses have also been observed in each of the 4 IRAC channels (3.6–8.0 $\mu$m) and the Infrared Spectrograph (IRS) at 16 $\mu$m (Stevenson et al. 2012). Knutson et al. (2009b) presented the first phase photometry, which spanned just over half the orbit in 8.0 $\mu$m. Most recently Zhang et al. (2018) published two full-phase observations in the Warm Spitzer bands (3.6 and 4.5 $\mu$m). We draw attention to §4.1 of Zhang et al. (2018), where the authors point out that inconsistencies between the bands warrant a fair degree of skepticism. In particular, the positive phase offset, or late minimum, of the 3.6 $\mu$m time series, is difficult to explain with modeling that assumes spin-orbit synchronization.

WASP-12 b was first discovered in transit via the SuperWASP camera (Hebb et al. 2009), and subsequent transits have been observed in the $V$ band (Chan et al. 2012), $R$ (Maciejewski et al. 2011), $J$, $H$, and $K_{s}$ (Croll et al. 2011), and at 3.6–8.0 $\mu$m from Spitzer (Campo et al. 2011; Cowan et al. 2012a; Stevenson et al. 2014). Mandell et al. (2013) also provided transit spectroscopy from WFC3 (1.1–1.7 $\mu$m) on the Hubble Space Telescope. Our work focuses on the full-orbit phase curves available from the warm Spitzer (3.6 and 4.5 $\mu$m) channels, originally published in Cowan et al. (2012a).

CoRoT-2 b has been observed in transit (Alonso et al. 2008) and eclipse via photometry in the Spitzer IRAC 3.6 (Deming et al. 2011), 4.5, and 8.0 $\mu$m (Gillon et al. 2010) channels. Transit spectra have also been measured (Bouchy et al. 2008; Czesla et al. 2012; Baluev et al. 2015). Refinement of the stellar parameters for CoRoT-2 led to revised transit depths (Southworth 2011, 2012; Baluev et al. 2015). Even with the revised radius, CoRoT-2 b is estimated to have a density comparable with Jupiter’s, hinting at a possible radius inflation that might be due to uncharacterized tidal processes. Recently, Dang et al. (2018) published a full-phase light curve at 4.5 $\mu$m, which showed an unusual and substantial $23^{\circ}$ westward phase offset. The authors proposed possible astrophysical sources of the offset including westward winds, magnetic effects, or partial cloud cover. Inhomogeneous cloud cover has been studied as a possible explanation of the westward offsets seen in optical light curves (Shporer 2017), and optical follow-up for CoRoT-2 b may support this hypothesis (Barstow 2018). Here we offer the alternative hypothesis of a high obliquity state. We show that this could not only produce the observed westward offset, but it would also provide new information about the formation and dynamical evolution mechanisms of close-in giants.

4 Components of the Oblique Thermal Model

In order to develop a model of thermal phase variations for a planet with non-zero spin obliquity, we build upon the model framework developed for the analysis in Adams & Laughlin (2018). Each planet is divided into a grid with cells of dimension $5^{\circ}\times 5^{\circ}$ by latitude and longitude. We start with the known system properties of each modeled planet (Table 2). There are 6 tunable parameters which jointly govern the resulting thermal emission. The first three — the albedo $A$, equilibrium radiative timescale $\tau_{\mathrm{rad}}$, and minimum temperature $T_{0}$ — most directly control the thermal properties of each cell. They are put into their formal context in §4.1. In short, each cell absorbs a fraction $1-A$ of radiation from the host star and re-radiates as a blackbody with a corresponding characteristic timescale. Its brightness temperature is set both by the time-dependent instellation and any non-stellar heating (e.g. tidal heat emanating from the planet interior), the latter of which is captured with the minimum temperature. The remaining three parameters define the components of the rotation vector. Movement of the cells in the model is set entirely by the rotation; the infrared photospheric layer of the planet rotates with some angular frequency $\omega_{\mathrm{rot}}$ ($\equiv 2\pi/P_{\mathrm{rot}}$ for the rotation period $P_{\mathrm{rot}}$) with some orientation of its axis relative to the orbital plane. This may be written in terms of two spherical angles: $\epsilon$, which is the angle between the planet’s axis and the orbit normal, and $\zeta$, which is the projected angle of the axis in the orbital plane relative to a reference direction. We set this reference direction along the planet-star line during periastron, such that $\zeta=0$ implies the northern hemisphere444The northern hemisphere is defined according to the axis around which the planet rotates counter-clockwise. Retrograde motion is therefore obtained by setting $\epsilon>90^{\circ}$. summer solstice occurs during periastron (equivalent to $f_{\mathrm{sol}}$ in Ohno & Zhang 2019a). This corresponds to an angle of the axial projection from the line of sight of $\zeta_{\mathrm{obs}}=\zeta-\nu_{\mathrm{tra}}$, for true anomaly at transit $\nu_{\mathrm{tra}}\equiv\pi/2-\varpi$; we assume $\varpi\rightarrow\pi/2$ for circular orbits so that $\zeta=\zeta_{\mathrm{obs}}$. The time evolution of grid cell temperatures is convolved with both the viewing geometry from Earth, which is approximated as edge-on, and the relevant instrumental band profiles, to predict full-orbit light curves for a given set of parameter values. Our fitting routine is detailed in §5.

4.1 Thermal Evolution of the Cells

The planetary model is initialized at apastron with a uniform surface temperature $T_{0}$; for planets on circular orbits, we set an arbitrary argument of periastron $\varpi=\pi/2$ such that apastron occurs during secondary eclipse. To calculate the incoming stellar radiation over the orbit, we start with the star-planet separation

[TABLE]

where $a$ is the orbital semi-major axis, $e$ the orbital eccentricity, and $\nu=\nu\!\left(t\right)$ the true anomaly. To solve for the true anomaly from the time in orbit we first calculate the mean anomaly $M$, which is directly proportional to time: $M\!\left(t\right)=\omega_{\mathrm{rot}}\left(t-t_{\mathrm{peri}}\right)$ for rotation rate $\omega_{\mathrm{rot}}$ and periastron passage time $t_{\mathrm{peri}}$. The eccentric anomaly $E=E\!\left(t\right)$ is then given by Kepler’s equation

[TABLE]

There is a direct relation between regular time intervals and regular intervals in mean anomaly, but not for the eccentric or true anomalies; we must calculate the latter numerically. The sine and cosine of the true anomaly, $\nu$ are given by

[TABLE]

At each time $t$ in the orbit and longitude/latitude $\left(\phi,\theta\right)$ on the planetary surface, the equilibrium temperature is calculated via

[TABLE]

where $A$ is the planetary albedo, $L_{\star}$ the stellar luminosity, $\sigma$ the Stefan-Boltzmann constant, $r=r\!\left(t\right)$ the star-planet separation, $\alpha_{\star}=\alpha_{\star}\!\left(\phi,\theta,t\right)$ the local stellar altitude, and $T_{0}$ the minimum temperature parameter. The stellar altitude $\alpha_{\star}$ is, relative to the sub-stellar point pointed to by $\hat{r}_{\star}=\phi_{\star}\hat{\phi}+\theta_{\star}\hat{\theta}$ and the unit normal $\hat{n}$ at the position,

[TABLE]

where $\hat{n}\cdot\hat{r}_{\star}\rightarrow\hat{r}\cdot\hat{r}_{\star}=\cos\theta\cos\theta_{\star}\left[\cos\left(\phi-\phi_{\star}\right)-1\right]+\cos\left(\theta-\theta_{\star}\right)$ for spherical planets, and the sub-stellar point is given by555For the equivalent derivation in inertial Cartesian coordinates, see §3 in Dobrovolskis (2013).

[TABLE]

For tidally distorted planets, the unit normal $\hat{n}$ no longer matches the position unit vector $\hat{r}$; §4.3 covers this case.

Finally, the change in temperature of each cell in time is calculated as

[TABLE]

which is a differential equation we can evaluate numerically for sufficiently small timesteps. We choose to divide each orbit into 200 timesteps.

4.2 Generating Observables

Given a temperature map $T\!\left(\phi,\theta,t\right)$, we solve for the corresponding planet-star flux contrast via

[TABLE]

where $B_{\lambda}\!\left(T\right)$ is the specific blackbody intensity at a wavelength $\lambda$ and temperature $T$, $w=w\!\left(\lambda\right)$ the weighted response of the instrumental bandpass at $\lambda$, and $V=V\!\left(\phi,\theta,t\right)$ is the component of the normal vectors of the cells along the line of sight, given by

[TABLE]

where the sub-observer point is given by

[TABLE]

and $\nu_{\mathrm{ecl}}=3\pi/2-\varpi$ is the true anomaly during secondary eclipse (i.e. at time $t_{\mathrm{ecl}}$), and for spherical planets, $\hat{n}\cdot\hat{r}_{\mathrm{obs}}=\cos\theta\cos\theta_{\mathrm{obs}}\left[\cos\left(\phi-\phi_{\mathrm{obs}}\right)-1\right]+\cos\left(\theta-\theta_{\mathrm{obs}}\right)$. However, as with the stellar altitude, for non-spherical planets the calculation is more involved, as we discuss in the following section.

4.3 A Simple Model of Tidal Distortion

A tidally distorted planet can have a significantly aspherical shape. We first outlined a model of tidal asphericity in Appendix A of Adams & Laughlin (2018). Here we adopt that work’s primary assumption to model the distorted shape as a prolate spheroid, with the long axis displaced clockwise in the orbital plane from the star-planet line by some lag angle $\lambda\equiv\cos^{-1}\!\left(\hat{r}_{\star}\cdot\hat{r}_{\ell}\right)$. We will assume this lag angle is zero for our analysis, but include it for completeness. The lengths of the long and short axes are dictated by the gravitational potential of the star-planet system. Consider a coordinate system where $\hat{z}$ points from the planet center along the long axis, $\hat{y}$ points along the orbit normal, and $\hat{x}$ points along the short axis according to a right-handed coordinate system. Then the potential becomes

[TABLE]

where $a$ is the orbital semi-major axis and $\xi$ is the planet-star mass ratio $M_{\mathrm{p}}/M_{\star}$. We then fit the cross-sectional area along $\hat{x}$ and $\hat{y}$ to the observed transit depth to get the long and short axis lengths666In general the solution yields unequal extents along $\hat{x}$ and $\hat{y}$, while the prolate spheroid assumption implies they are equal. For WASP-12 b, the differences are minor enough that the corresponding differences in calculated ellipticity are $\lesssim$ 0.01..

For an ellipsoid the area of each cell is position-dependent, and will therefore affect the surface area over which it radiates. Once we have the planetary semi-major and semi-minor extents (defined as one-half of the long and short axes, respectively), denoted $A_{\mathrm{p}}$ and $B_{\mathrm{p}}$, we can calculate the areas of individual cells. We adapt the result from Equation A7 of Adams & Laughlin (2018) for the area of a cell spanning longitudes $\phi_{i}$–$\phi_{j}$ and latitudes $\theta_{i}$–$\theta_{j}$, now with a more complicated relationship between the Cartesian coordinates for the planet:

[TABLE]

where $\chi\equiv A_{\mathrm{p}}/B_{\mathrm{p}}$ is the axis ratio and

[TABLE]

A second effect of the non-spherical shape is a change in the stellar altitude as a function longitude and latitude. To quantify the change we first write the function representing the shape of the prolate ellipsoid:

[TABLE]

We need a way of expressing these ellipsoidal coordinates in the oblique coordinates (i.e. with respect to the latitude/longitude coordinates defined by the rotation). To do this we note that we can express the oblique positions of both the sub-stellar point, given by equation 6, and the extreme point of the planet along the long axis, given by

[TABLE]

If the lag angle $\lambda>0$, then we can construct our ellipsoidal Cartesian unit vectors entirely with respect to the unit position vectors for these two points:

[TABLE]

From this the normal unit vector at a given point on the surface is given by

[TABLE]

where

[TABLE]

Then the cosines of the stellar altitude and the visibility are given by the component of $\hat{n}$ along the instellation and observer lines, as in Equations 5 and 9.

[TABLE]

and the component along the observer line is

[TABLE]

5 Fit Methods and Results

For each planet we run two cases: one with 6 free parameters where the rotation rate is unconstrained, and another with 5 parameters where the rotation rate is fixed to the obliquity-dependent equilibrium rotation rate

[TABLE]

where $n\equiv 2\pi/P_{\mathrm{orb}}$ is the mean motion (see §7.2.1 for more details). Additionally, for HD 149026 b and WASP-12 b, which have data in two bands, we run the above two cases where we constrain the axial orientation to be consistent between the bands, for a total of 4 oblique model cases.

We use a Markov-Chain Monte Carlo process to broadly evaluate the likelihood landscape over the relevant parameters and converge on the sets of parameter values for each case with the most favorable likelihoods. Specifically, we employ a Metropolis-Hastings algorithm with simulated annealing (Kirkpatrick et al. 1983). Annealing introduces a temperature parameter that corresponds to acceptance probability, and serves to both broadly explore the likelihood space at initially high values, and to converge on optimal solutions as it is gradually reduced. After convergence, we continue the MCMC chain without annealing, to estimate uncertainties on the parameter values. The quoted 1-$\sigma$ uncertainties are determined by the 68% ranges on either side of the values. The best-fit values and uncertainties for the models with band-distinct obliquities are listed in Table 2. The resulting model light curves are plotted with these 1-$\sigma$ uncertainties, as well as the 2-$\sigma$ uncertainties, determined by the 95% ranges, in the leftmost columns of Figures 2–4.

To evaluate the relative fit quality of these oblique models, we calculate the Akaike Information Criteria (Akaike 1973, 1974) with a second-order correction for small sample sizes (Sugiura 1978; M. Hurvich & Tsai 1989, 1995). This is given by

[TABLE]

where $k$ is the number of free parameters, $\ln\mathcal{L}$ the log likelihood (and, by extension, the full second term $-2\ln\mathcal{L}$ being equal to the chi-squared statistic), and $n$ the sample size. We express the AIC values relative to the smallest value (i.e. the value of the most favorable model) in each band ($\Delta\mathrm{AICc}\equiv\mathrm{AICc}-\mathrm{AICc}_{\mathrm{min}}$), in Table 3. Additionally, to more appropriately compare the models with band-distinct obliquities with the models co-varying in obliquity, we calculate the cumulative AICc values for the combined model across both bands. Following the interpretation of Burnham & Anderson (2004), $\Delta\mathrm{AICc}\leq 2$ indicates a substantial level of evidence for a model (as compared with the most favorable); $4\leq\Delta\mathrm{AICc}\leq 7$ implies low evidence, and $\Delta\mathrm{AICc}>10$ effectively implies no evidence. We discuss the results of each planet in the following sub-sections.

5.1 HD 149026 b

Following the prescription of the AIC, we find that the preferred model at 3.6 $\mu$m is the model with zero obliquity and unconstrained rotation. This implies that any fit improvement from adding axial orientation parameters are not statistically warranted. Indeed, the least-constrained model, which has a distinct obliquity and unconstrained rotation, returns a nearly upright axis, effectively approximating the non-oblique case. The $\Delta\mathrm{AIC}\sim 4$ is consistent with nearly identical likelihoods, since we are increasing the number of parameters by 2 by including obliquity.

In contrast, the most preferred model at 4.5 $\mu$m is the simplest — zero obliquity and synchronous rotation. This implies that any discernible phase offset in the data is not strong enough to warrant a model that can capture it, either by a sub-synchronous rotation or high obliquity. The non-oblique, unconstrained rotation case is next in line, with a best-fit spin-orbit ratio of nearly two. Moving into the oblique cases, we see that there is a moderate drop-off in support for the models at 4.5 $\mu$m; this is not surprising, since the light curves at 3.6 and 4.5 $\mu$m disagree with the direction of offset.

Interestingly, when we compare the combined model criteria, the combined non-oblique model with unconstrained rotation is on top, but the model with consistent obliquity and sub-synchronous rotation is not far behind. The most likely interpretation is that the 3.6 $\mu$m data are modestly better fit by an axis nearly perpendicular to the orbit normal, rather than simply a slow effective rotation rate. Given the very slow returned rotation for this case, it also suggests that the rotation has a weak effect on the offset at this high of obliquity, since the transverse motion is taken almost completely out of the plane of the sky from our perspective.

5.2 WASP-12 b

Due to its strong tidal distortion, WASP-12 b is the most complicated of the planets, and we expect that it will be difficult to fully capture its shape only accounting for rotation and obliquity. In sharp contrast with the moderate differences in $\Delta\mathrm{AIC}$ between models in the bands of HD 149026 b, which imply a dominant effect of the number of parameters, for WASP-12 b both the synchronous and sub-synchronous models are vastly worse than those with unconstrained rotation. This makes sense given that the light curves in both bands show strong variations between occultations, at least one of which can be approximated using an eastward phase offset. Eastward offsets are more amenable to a super-rotating atmospheric layer than high obliquity, and indeed the models with the most free parameters return small values for obliquity.

AIC however cannot provide a quantification of the absolute fit quality. While our model light curves generally lie within the range of observed fluxes at any given time, the maximum observed amplitude of semi-annual variations exceed what can be captured with a combination of spin-orbit geometry and thermal modeling. Ellipsoidal variations due to the planet’s tidal distortion should only change the observed surface area by $\sim\!10\%$ (Li et al. 2010; Lai et al. 2010), and the complementary effects on the star’s shape should be a further order of magnitude smaller. These correspond to variations on the orders of $\sim\!4\times 10^{-4}$ and $\sim\!4\times 10^{-5}$, according to Cowan et al. (2012a). Indeed, the incorporation of an ellipsoidal distortion to the non-oblique thermal model (e.g. in Adams & Laughlin 2018) does not provide enough additional variation to approach the full amplitudes of the data.

5.3 CoRoT-2 b

Our analysis of the light curves is more limited for CoRoT-2 b since it only has one band, but the AIC values do indicate that the non-synchronous models are much better fits to the data. As with the 3.6 $\mu$m fits for HD 149026 b, the differences in selection criteria between the non-oblique, non-synchronous and oblique models appear to be largely influenced by the AIC’s penalty for additional parameters by introducing axial orientation. From this we infer that the model fits from either a westward wind interpretation (slower than synchronous rotation with little to no obliquity) or a sub-synchronous rotation at moderately high obliquity are of similarly effective quality.

Considering these three cases, we can draw some conclusions. For planets whose phase variations are already amenable to fitting by a non-oblique thermal model, but whose observed westward phase offsets require a slower-than-synchronous effective rotation rate (most readily interpreted as westward winds), our fit quality is at least as good with an oblique model where we assume an obliquity-dependent sub-synchronous rotation (Equation 21). However, particularly for WASP-12 b, we have no cases where the addition of obliquity improves the quality of fit of a thermal model which has difficulty reproducing all major features of the data.

6 Dynamical Evolution of a Tightly Packed System

Having established the feasibility of obliquity to fit the observed phase variations, we now place this interpretation in a dynamical framework in the following sections. Batygin et al. (2016) set a physical framework where a system comprising multiple close-in planets (i.e. $P\lesssim 100$ days) with masses in the super-Earth regime ($1<M/M_{\oplus}<30$) can evolve to have high mutual orbital inclinations. They first determine that, above an initial mass of $\gtrsim 15M_{\oplus}$, the innermost planet can undergo runaway accretion under a range of stellar nebular densities, reaching $M_{p}\sim M_{\mathrm{J}}$ within a timescale $\sim\!10^{6}$ years, thereby predicting in-situ formation of close-in giant planets. Their model couples this with the evolution of the host star onto the main sequence and determines that outer planets can evolve to high mutual inclinations via nodal regression commensurability. These results constitute a prediction that systems with hot Jupiters and mutually inclined outer planets should be common.

In order to evaluate the feasibility of evolving a dynamically full multi-planet system into a configuration consistent with the existence of secular spin-orbit resonance and associated high obliquities, we performed dynamical simulations using the REBOUND777https://github.com/hannorein/rebound software package with the IAS15 integrator (Rein & Liu 2012; Rein & Spiegel 2015). For our fiducial precursor system, we adopted a known multi-planet system, Kepler-107 (Rowe et al. 2014; Van Eylen & Albrecht 2015). We then statistically explored the dynamical reaction of the configuration to orbital instability triggered by rapid increase in the mass of the innermost planet. We assumed an instantaneous increase from the estimated innermost planet mass of 3.7 $M_{\oplus}$, to that of Jupiter. We also integrated the system without any mass increase. For both cases we performed 100 separate simulations with a duration of at least $\sim\!10^{5}$ Earth years ($\sim\!10^{7}$ orbits of the innermost planet as initialized).

The orbital geometry of the Kepler-107 system is not entirely constrained by the transit parameters reported by Rowe et al. (2014). We therefore adopted initial values for the semi-major axes from Rowe et al. (2014) and initial eccentricities from Van Eylen & Albrecht (2015). For the angular orbital elements, we used the known transit parameters in conjunction with the following process to randomly generate orbital elements consistent with the observational constraints.

Set the line of sight to the positive $x$-axis, the default reference direction in REBOUND888Note that in this setup, an inclination of zero will imply the orbit is viewed edge-on, rather than the standard $i=90^{\circ}$. We adopt this non-standard definition purely for convenience.. 2. 2.

Choose a longitude of periastron, $\varpi\sim U\!\left(0,2\pi\right)$. 3. 3.

Choose a longitude of ascending node $\Omega$ from a uniform distribution over the set of angles satisfying

[TABLE]

where $b$ is the impact parameter, $R_{\star}$ and $R_{\mathrm{P}}$ the stellar and planetary radii, respectively, and

[TABLE]

the star-planet separation at the transit mid-point. 4. 4.

Calculate the inclination in the range $\left[0,\pi/2\right]$ that solves the transit condition

[TABLE]

After integrating systems constructed in this manner for $\sim\!10^{5}$ Earth years, we see a significant evolution of the systems away from their initially nearly co-planar states. We calculate the inclinations (relative to the net system angular momentum vector) for the final orbits of the remaining planets and evaluate the transit conditions along the initial line of sight (Figure 5). In a significant fraction of simulations, some combination of the outer planets ($c$–$e$) no longer transit due to large inclinations which can exceed 20∘. Out of the 100 simulations with the nominal Kepler-107 masses, only 3 saw an ejection (defined in our simulation is defined as reaching a separation from the star an order of magnitude larger than the initial semi-major axis). The fraction of ejections increases greatly in the inner-Jupiter model, with a significant fraction of trials ejecting $c$ and/or $d$ within $10^{5}$ years. Half of ejections occurred within the first $\sim\!2.5\times 10^{4}$ years. These results as a whole suggest that, under a significant perturbation such the runaway accretion of the inner planet, it is plausible that a compact multi-transiting planet system such as Kepler-107 could evolve to high mutual inclinations and/or ejections for the outer planets, leaving a majority of cases where only the inner planet could be observed to transit.

The simulations presented here have only accounted for the radii of the star and orbiting planets in terms of collision detection, whose prescribed resolution was a merger. No mergers occurred in any trial. Otherwise, the bodies are effectively treated as point masses. If these bodies are endowed with structure, as formulated for example by Mardling & Lin (2002), then secular inclination resonances can act to increase the degree of dynamical instability (leading to agglomerating collisions) and mutual inclinations among remaining planets, thereby increasing the fractions of non-transiting planets (Batygin et al. 2016).

7 Cassini State Driven by an Inclined, External Perturber

While we have examined the possibility that close-in giant planets have non-zero obliquities, we have not yet discussed in detail how such states may arise. Tidal torques are strongest for these close-in, large-radius planets, and they act to dampen planetary obliquities to zero. Large tilts can be maintained, however, if a planet is locked in a secular spin-orbit resonance involving synchronous precession of its spin vector and orbital angular momentum vector (Fabrycky et al. 2007). Of particular interest is Cassini State 2 (Peale 1969), where the planet maintains a large obliquity as the spin and orbital axes precess at the same rate on opposite sides of the axis normal to the invariable plane. In the absence of dissipation, these three axes are coplanar, but in the presence of tides, the spin axis is slightly shifted out of the plane.

In this section, we evaluate the possibility that the planets HD 149026 b, WASP-12 b, and CoRoT-2 b have their spin axes trapped in Cassini states due to spin-orbit resonances driven by exterior perturbing planets. The case of a resonance for WASP-12 b was already discussed in detail by Millholland & Laughlin (2018) as a theory for the planet’s rapid orbital decay (Maciejewski et al. 2016; Patra et al. 2017; Maciejewski et al. 2018; Bailey & Goodman 2019). We include elements of Millholland & Laughlin (2018)’s analysis here for consistency and comparison with the other two systems.

We start by calculating the frequencies of spin and orbital precession. We then examine the constraints on the hypothetical perturbing planets in order for the spin-orbit resonant configurations to be plausible. Table 7.2.1 shows the system parameters we adopt in these calculations.

7.1 Spin-Orbit Resonant Frequencies

The torque from the host star on a rotationally-flattened planet will cause the planet’s spin-axis to precess about the orbit normal at a period, $T_{\alpha}=2\pi/(\alpha\cos\epsilon)$. Here $\epsilon$ is the obliquity and $\alpha$ is the precession constant (Ward & Hamilton 2004), which is given by

[TABLE]

We have assumed there are no satellites and have defined $M_{\mathrm{p}}$, the planet mass, $C$, the moment of inertia normalized by $M_{\mathrm{p}}{R_{\mathrm{p}}}^{2}$, and $\omega$, the spin frequency. This expression also assumes that the coefficient of the quadrupole moment of the planet’s gravitational field, $J_{2}$, takes the form (Ragozzine & Wolf 2009),

[TABLE]

The secular spin-orbit resonance requires a commensurability between the planet’s spin-axis precession frequency and its orbit nodal regression frequency, $g=\dot{\Omega}$. Planet-planet interactions are one source of nodal regression. In a two-planet system, the nodes of both planets regress uniformly due to secular perturbations. The frequency is given by Laplace-Lagrange theory to be

[TABLE]

if the planets are not near mean-motion resonance (Murray & Dermott 1999). Here, $\alpha_{12}=a_{1}/a_{2}$ and $n_{i}$ is the mean-motion of planet $i$, $n_{i}^{2}=GM_{\star}/a_{i}^{3}$. The constant, $b_{3/2}^{(1)}(\alpha_{12})$ is a Laplace coefficient, defined by

[TABLE]

In addition to planet-planet interactions, the stellar quadrupole gravitational moment also induces orbit nodal recession about the stellar spin vector. In the absence of secular planet interactions, this occurs at the frequency (Spalding & Batygin 2017)

[TABLE]

Accordingly, in a multiple-planet system, the planets’ orbit normal vectors evolve in response to perturbation components at different frequencies. One frequency, $g_{p-p}$, is due to the planet-planet perturbations and the other, $g_{\star}$, is associated with the stellar quadrupole moment. These are close but not exactly equal to the equations 28 and 30 above, since the analytical expressions only account for one driver of nodal recession at a time. When the stellar equatorial plane is coincident with the plane perpendicular to the total orbital angular momentum vector, these two frequencies add linearly such that $g=g_{p-p}+g_{\star}$, and the nodal recession is uniform. However, if these planes are not coincident (i.e. if there is an angle between the stellar spin vector and the total orbital angular momentum vector), then the frequencies do not add linearly, but rather the nodal recession is a non-uniform superposition of these modes.

The spin-orbit resonance can be encountered and captured when the spin precession is commensurable with either one of these orbital frequency components, $g_{p-p}$ or $g_{\star}$ (Millholland & Laughlin 2019). In the case of short-period planets, however, the spin axis precession is so fast that spin-orbit resonances induced by planet-planet interactions are much more likely. Although $g_{\star}$ is fast early in the system’s lifetime because the star is rapidly rotating and has not finished contracting (Batygin & Adams 2013; Spalding & Batygin 2017), $g_{\star}$ is significantly smaller for a main-sequence star, and it is generally not fast enough for spin-orbit resonance. To illustrate this explicitly, we calculate $g_{\star}$ for our three case studies. The results are shown in Table 7.2.1. In all three systems, $g_{\star}$ is much smaller than the spin axis precession constant. This clearly indicates that nodal recession must be provided by planet-planet perturbations if these systems are in spin-orbit resonances. We will therefore assume in the remainder of this work that planet-planet resonances are the relevant ones.

7.2 Constraints on the masses and semi-major axes of potential inclined companions

Suppose that, in each of the three systems, planet $b$ is in a high-obliquity spin-orbit resonance with its orbital precession induced by secular interactions with an exterior, inclined, and as-yet undetected companion. Existence of this configuration places significant constraints on the masses, semi-major axes, and inclinations of the perturbing companions. These constraints stem from

(i) upholding the Cassini state, (ii) preserving total angular momentum conservation, and (iii) maintaining consistency with existing radial velocity (RV) data. In the subsections that follow, we review and apply each of these constraints in detail.

7.2.1 Constraints from Maintaining the Cassini State

If a planet is captured in a Cassini state, there will be a resonant commensurability that can be stated as (Ward & Hamilton 2004)

[TABLE]

This holds if the planetary obliquity, $\epsilon$, is large compared to the mutual inclination between the orbits. We will take that as an assumption and use this condition to calculate the range of values of each perturbing planet’s mass, $M_{p2}$, and semi-major axis, $a_{2}$, that allow for resonant commensurability. Table 7.2.1 shows the estimates of the spin axis precession periods used in these calculations.

In Figure 6, we show heatmaps in $M_{p2}$ and $a_{2}$ space that represent the obliquity of planet $b$ necessary for the resonance to hold in each system. We assume that the nodal recession is driven by interactions with the secondary planet, such that $g=g_{p-p}\approx g_{{}_{\mathrm{LL}}}$. We also assume $e_{b}=0$ and that the spin rate of planet $b$ is at equilibrium, at which $d\omega/dt=0$. The equilibrium rate is given by equation 21 in the traditional viscous approach to equilibrium tide theory (Levrard et al. 2007). The quantity $n\equiv 2\pi/P_{\mathrm{orb}}$ is the mean-motion. Combining equations 26, 31, and 21, the resonance condition becomes

[TABLE]

where $\alpha_{\mathrm{syn}}$ is the value of $\alpha$ in the case of synchronous rotation, $\omega=n$. The solution for $\epsilon$ is then

[TABLE]

This expression can be used to calculate $\epsilon$ for a range of values of $M_{p2}$ and $a_{2}$.

Figure 6 shows that if $a_{2}$ is too small, $\lvert g|$ is too large and no resonance is possible (white regions). Alternatively, if $a_{2}$ is too large, $\lvert g|$ is small and $\epsilon\sim 90^{\circ}$, which is unstable in the long-term. Figure 6 thus allows us to define approximate constraints on the semi-major axes of the hypothetical companion planets. The constraints on their masses are not strong due to the weak dependence of $g_{{}_{\mathrm{LL}}}$ on $M_{p2}$ when $M_{p1}\gg M_{p2}$.

7.2.2 Constraints from Angular Momentum Conservation

For the equilibrium tidal theory that we are considering, a planet maintained in an oblique state spins sub-synchronously (equation 21), and the rate at which tides convert orbital energy into heat energy is orders of magnitude larger than it would be in the case of zero obliquity (Levrard et al. 2007; Wisdom 2008). An oblique tidally dissipating planet migrates inwards towards the star, decreasing $\lvert g\rvert/\alpha$ in the process and further increasing its obliquity. The decrease in orbital angular momentum associated with inward migration must be counteracted such that total angular momentum of the system, $\mathbf{J}$, is conserved. This can be accomplished through gradual alignment of the planetary orbital angular momenta, $\mathbf{L_{1}}$ and $\mathbf{L_{2}}$, with one another and with the stellar spin angular momentum, $\mathbf{S_{\star}}$. If we assume for a conservative argument that alignment with $\mathbf{S_{\star}}$ does not play a role, then the second planet must have enough angular momentum to preserve the total conservation. This allows conservative limits to be placed on the orbits of the perturbing planets (Fabrycky et al. 2007).

We begin with an expression for $\mathbf{J}$. We define $i$ to be the angle between $\mathbf{L_{1}}$ and $\mathbf{L_{2}}$. In addition, we define $\phi$ to be the angle between $\mathbf{S_{\star}}$ and $\mathbf{L_{1}}+\mathbf{L_{2}}$. Then, assuming the planets’ own spin angular momenta are negligible, the magnitude of $\mathbf{J}$ is given by

[TABLE]

Secular interactions between the planets do not change $a_{2}$ or $e_{2}$ to first order, so $L_{2}$ remains fixed. Therefore, as $L_{1}$ decreases, conservation of $J$ must be upheld by an increase in $S_{\star}$ via tidal spin-up of the star (Brown et al. 2011) or a decrease in $i$ and $\phi$ through reorientation of the orbits. As stated above, we will conservatively assume that stellar spin-up and reorientation between the orbital and stellar spin angular momenta do not play a role. Under this assumption, the maximum possible value of planet $b$’s initial semi-major axis, $\max({a_{1i}})$, may be expressed in terms of $a_{2}$ and $M_{p2}$,

[TABLE]

where $a_{1}$ corresponds to the present-day value of the semi-major axis. This expression is obtained by assuming that $i$ was initially near $90^{\circ}$ and is currently near $0^{\circ}$. We also assumed $e_{2}\approx 0$.

This expression does not hold enough information to constrain $a_{2}$ and $M_{p2}$ in and of itself, since there is no clear limit on $\max(a_{1})$. To develop a constraint, we apply the additional requirement that the initial obliquity of planet $b$ must be near $0^{\circ}$. Though this is not strictly necessary, if it is true, it makes the initial resonant capture scenario easy to explain. (Recall that $\epsilon$ increases as the planet tidally migrates inwards.) For a given $a_{2}$ and $M_{p2}$, the semi-major axis, $a_{1i}$, at which the initial obliquity is zero satisfies the expression

[TABLE]

There is only a plausible solution if $a_{1i}<=\max({a_{1i}})$. Accordingly, this is the additional constraint that we use in conjunction with equation 35.

The green lines in Figure 6 delineate the region that simultaneously upholds the conservative angular momentum limits and maintains the possibility of a small initial obliquity for planet $b$. Significant regions of phase space are ruled out for both the WASP-12 and CoRoT-2 systems. While this is somewhat problematic for the theory, it is important to keep in mind that these tight constraints can be strongly alleviated and/or removed by assuming a non-zero initial obliquity or by allowing realignment with $\mathbf{S_{\star}}$ to account for some degree of the system’s angular momentum conservation.

7.2.3 Constraints from Radial Velocity Data

An entirely unrelated constraint arises from the set of RV measurements of the systems. The hypothetical companion planets must have RV semi-amplitudes smaller than the current detection limits; otherwise they would have already been discovered. Table 7.2.3 shows current estimates of the RV semi-amplitudes, $K_{1}$, of the known planets and the RV jitters of the fits, $\sigma_{\mathrm{jit}}$. The hypothetical outer planets must have $K_{2}\lesssim\sigma_{\mathrm{jit}}$. Figure 6 shows that the parameter space in the HD 149026 and WASP-12 systems are somewhat constrained by the RV limits. CoRoT-2’s jitter, however, is so large that hardly any phase space area is ruled out.

Given its small $\sigma_{\mathrm{jit}}$, it is worthwhile to examine the RVs of HD 149026 to search for any hints of a second signal. The Lick-Carnegie Exoplanet Survey Team (LCES) collected 70 measurements of Keck/HIRES Doppler velocities (Butler et al. 2017) over 8.5 years. The top panel of Figure 7 shows the Keck RVs phase-folded at the period of planet $b$, $P=2.8759$ days. The second panel shows the residuals obtained after fitting a circular orbit. Finally, the bottom panel displays a Lomb-Scargle periodogram of the RV residuals. There is a noticeable peak at 12.68 days. This peak may be related to stellar activity, since it is close to the $\sim\!13$ day stellar rotation period suggested by the $v\sin i$ measurement from Sato et al. (2005). Alternatively, it may be the signature of an additional planet, and a closer examination of both the rotational signature in the existing spectra, as well as additional Doppler measurements of the star may be warranted.

7.3 Summary of Constraints

In summary, the combination of the requirements of a secular spin-orbit resonance for planet b, total angular momentum conservation, and RV detection limits place strong constraints on the parameters of hypothetical perturbing planets. The limits from angular momentum conservation are particularly restrictive, though we used very conservative assumptions in those calculations. All three systems could therefore host additional, as-yet-undetectable planets with parameters appropriate for generating high-obliquity Cassini states for their companion giant planets.

8 Discussion

While most thermal full-phase light curves of close-in giant planets show that the hottest region on the planet is eastward of the sub-stellar point, which is consistent with super-rotating winds, there are now at least two planets (HD 149026 b and CoRoT-2 b) with significant westward hotspot offsets in at least one Spitzer band (Zhang et al. 2018; Dang et al. 2018). One way of reconciling these westward offsets is to consider variations in the spin axis orientation, which fundamentally changes the relationship between the instellation variations due to the rotation and those due to the orbit. Planets on very short-period orbits occupy an interesting regime where we expect the rotation and orbital rates to be comparable; changes in the spin geometry can therefore have major effects on the observed phases. In this work, we developed a thermal radiative model to investigate how non-zero planetary obliquities may produce observable signatures in full-orbit phase curves.

It is important to note that we have not considered the degeneracies inherent to a modeling approach that considers only the output 1-D photometry. Many advances have been made recently that focus on first constructing physical 2- or 3-dimensional maps of the planetary flux, then calculating photometric variations (e.g.  Rauscher et al. 2018). These will likely be the most rigorous frameworks for future photometric modeling. Despite these limitations, and acknowledging that there are multiple potential explanations for westward offsets, we showed that, to first order, highly-oblique hot Jupiters should have different phase variation morphologies from those with little to no obliquity.

We focused on a case study of three planets — HD 149026 b, WASP-12 b, and CoRoT-2 b — all of which have previous light-curve observations that exhibit anomalous features (Cowan et al. 2012a; Zhang et al. 2018; Dang et al. 2018; Barstow 2018) that we propose may arise from non-zero obliquities. We reanalyzed the existing Spitzer photometry using a thermal model that allows the spin orientation to vary freely.

A combination of sufficiently long thermal timescales and spin axes nearly perpendicular to the orbit normal can reproduce the strong westward offset in the 3.6 $\mu$m data of HD 149026 b. However, these oblique models are not preferred over a non-oblique model with fewer parameters, where the westward offset comes from an effective rotation slower than synchronous. The results are quite similar for CoRoT-2 b, where offset-capable models are strongly preferred, and oblique models can effectively reproduce the offset with a very tight constraint on the axial orientation, but due to the slight increase in the number of parameters are statistically not preferred over the simpler non-oblique model. For WASP-12 b, variations near quadrature cannot be effectively reproduced by a simple thermal model even with free spin geometry, though they fare better than models with constrained rotation. This suggests that more complex models with tidal distortion or refined instrumental characterization need to be considered to help explain the phase curves of WASP-12 b.

Among the three planets we studied, CoRoT-2 b has the tightest constraints on its proposed obliquity, at $45.8^{\circ}\pm 1.4^{\circ}$. In addition to the westward hotspot offset (Dang et al. 2018), a high obliquity state could also account for the planet’s extremely inflated radius, which is more anomalous than typical hot Jupiters (Guillot & Havel 2011). The enhanced tidal dissipation that accompanies a high obliquity state is more than enough to provide the heat. To second order in eccentricity, the rate of tidal dissipation in a state of equilibrium rotation and according to equilibrium tide theory is (Levrard et al. 2007)

[TABLE]

Here $K$ is given by

[TABLE]

Using $k_{2}=0.1$ and $C=0.2$ as in section 7, $Q=5\times 10^{5}$, $e=0$, and $\epsilon=45.8^{\circ}$ as suggested by the sub-synchronous fit, the tidal dissipation rate is $\dot{E}_{\mathrm{tide}}=1.2\times 10^{28}\ \mathrm{erg\ s^{-1}}$. This rate is a substantial fraction of the $2.7\times 10^{29}\ \mathrm{erg\ s^{-1}}$ of energy from incident stellar radiation, which indicates that obliquity tides may provide a clean solution to the anomalous radius inflation. Although the thermal phase curve modeling of WASP-12 b was inconclusive due to tidal distortion complications, this planet also has a highly inflated radius, which might similarly be a product of obliquity tides (Millholland & Laughlin 2018). Variations in the measured optical eclipse depths for WASP-12 b point to a possible additional occultation from escaping material (von Essen et al. 2019), which could correspondingly cause unexpected thermal variability.

If any of these planets do indeed have non-zero obliquities, these configurations could be induced by secular spin-orbit resonances with as-yet undetected exterior planets. We outlined dynamical arguments for each system to identify the regions of parameter space where such obliquity-producing perturbing planets could exist. There is allowable parameter space in each system, particularly in the HD 149026 system. Such results leave open the possibility that these close-in planets may be locked in a stable, high-obliquity spin states.

S.M. is supported by the National Science Foundation Graduate Research Fellowship Program under Grant Number DGE-1122492. This material is based upon work supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under Cooperative Agreement Notice NNH13ZDA017C issued through the Science Mission Directorate. We acknowledge support from the NASA Astrobiology Institute through a cooperative agreement between NASA Ames Research Center and Yale University. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Exoplanet Archive, Kepler, Spitzer (IRAC)