# Spi-OPS: Spitzer and CHEOPS confirm the near-polar orbit of MASCARA-1 b and reveal a hint of dayside reflection

M. J. Hooton, S. Hoyer, D. Kitzmann, B. M. Morris, A. M. S. Smith, A. Collier Cameron, D. Futyan, P. F. L. Maxted, D. Queloz, B. -O. Demory, et al.

► To cite this version:

M. J. Hooton, S. Hoyer, D. Kitzmann, B. M. Morris, A. M. S. Smith, et al., Spi-OPS: Spitzer and CHEOPS confirm the near-polar orbit of MASCARA-1 b and reveal a hint of dayside reflection. Astronomy and Astrophysics - A&A, 2022, 658, 10.1051/0004-6361/202141645. insu-03647282

HAL Id: insu-03647282

https://insu.hal.science/insu-03647282

Submitted on 20 Apr 2022

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
# Spi-OPS: *Spitzer* and CHEOPS confirm the near-polar orbit of MASCARA-1 b and reveal a hint of dayside reflection*

M. J. Hooton¹, S. Hoyer², D. Kitzmann³, B. M. Morris³, A. M. S. Smith⁴, A. Collier Cameron⁵, D. Futyan⁶,
P. F. L. Maxted⁷, D. Queloz⁶,⁸, B.-O. Demory³, K. Heng³,⁹, M. Lendl⁶, J. Cabrera⁴, Sz. Csizmadia⁴,
A. Deline⁶, H. Parviainen¹⁰, S. Salmon⁶, S. Sulis², T. G. Wilson⁵, A. Bonfanti¹¹, A. Brandeker¹²,
O. D. S. Demangeon¹³,¹⁴, M. Oshagh¹⁰, C. M. Persson¹⁵, G. Scandariato¹⁶, Y. Alibert¹, R. Alonso¹⁰,¹⁷,
G. Anglada Escudé¹⁸,¹⁹, T. Bárczy²⁰, D. Barrado²¹, S. C. C. Barros¹³,¹⁴, W. Baumjohann¹¹, M. Beck⁶,
T. Beck¹, W. Benz¹,³, N. Billot⁶, X. Bonfils²², V. Bourrier⁶, C. Broeg¹,³, M.-D. Busch¹, S. Charnoz²³,
M. B. Davies²⁴, M. Deleuil², L. Delrez²⁵,²⁶,⁶, D. Ehrenreich⁶, A. Erikson⁴, J. Farinato²⁷, A. Fortier¹,³, L. Fossati¹¹,
M. Fridlund²⁸,²⁹, D. Gandolfi³⁰, M. Gillon²⁵, M. Güdel³¹, K. G. Isaak³², K. Jones³, L. Kiss³³,³⁴,³⁵,
J. Laskar³⁶, A. Lecavelier des Etangs³⁷, C. Lovis⁶, A. Luntzer³⁸, D. Magrin³⁹, V. Nascimbeni³⁹, G. Olofsson¹²,
R. Ottensamer³⁸, I. Pagano¹⁶, E. Pallé¹⁰, G. Peter⁴⁰, G. Piotto³⁹,⁴¹, D. Pollacco⁹, R. Ragazzoni³⁹,⁴¹, N. Rando⁴²,
F. Ratti⁴², H. Rauer⁴,⁴³,⁴⁴, I. Ribas¹⁸,¹⁹, N. C. Santos¹³,¹⁴, D. Ségransan⁶, A. E. Simon¹, S. G. Sousa¹³,
M. Steller¹¹, Gy. M. Szabó⁴⁵,⁴⁶, N. Thomas¹, S. Udry⁶, B. Ulmer⁴⁷, V. Van Grootel²⁶, and N. A. Walton⁴⁸

*(Affiliations can be found after the references)*

Received 26 June 2021 / Accepted 10 September 2021

## ABSTRACT

**Context.** The light curves of tidally locked hot Jupiters transiting fast-rotating, early-type stars are a rich source of information about both the planet and star, with full-phase coverage enabling a detailed atmospheric characterisation of the planet. Although it is possible to determine the true spin-orbit angle $\Psi$ – a notoriously difficult parameter to measure – from any transit asymmetry resulting from gravity darkening induced by the stellar rotation, the correlations that exist between the transit parameters have led to large disagreements in published values of $\Psi$ for some systems.

**Aims.** We aimed to study these phenomena in the light curves of the ultra-hot Jupiter MASCARA-1 b, which is characteristically similar to well-studied contemporaries such as KELT-9 b and WASP-33 b.

**Methods.** We obtained optical Characterising ExOPlanet Satellite (CHEOPS) transit and occultation light curves of MASCARA-1 b, and analysed them jointly with a *Spitzer/IRAC* 4.5 µm full-phase curve to model the asymmetric transits, occultations, and phase-dependent flux modulation. For the latter, we employed a novel physics-driven approach to jointly fit the phase modulation by generating a single 2D temperature map and integrating it over the two bandpasses as a function of phase to account for the differing planet-star flux contrasts. The reflected light component was modelled using the general ab initio solution for a semi-infinite atmosphere.

**Results.** When fitting the CHEOPS and *Spitzer* transits together, the degeneracies are greatly diminished and return results consistent with previously published Doppler tomography. Placing priors informed by the tomography achieves even better precision, allowing a determination of $\Psi = 72.1^{+2.5}_{-2.4}$ deg. From the occultations and phase variations, we derived dayside and nightside temperatures of $3062^{+66}_{-68}$ K and $1720 \pm 330$ K, respectively. Our retrieval suggests that the dayside emission spectrum closely follows that of a black-body. As the CHEOPS occultation is too deep to be attributed to blackbody flux alone, we could separately derive geometric albedo $A_g = 0.171^{+0.066}_{-0.068}$ and spherical albedo $A_s = 0.266^{+0.097}_{-0.100}$ from the CHEOPS data, and Bond albedo $A_B = 0.057^{+0.083}_{-0.101}$ from the *Spitzer* phase curve. Although small, the $A_g$ and $A_s$ indicate that MASCARA-1 b is more reflective than most other ultra-hot Jupiters, where H$^-$ absorption is expected to dominate.

**Conclusions.** Where possible, priors informed by Doppler tomography should be used when fitting transits of fast-rotating stars, though multi-colour photometry may also unlock an accurate measurement of $\Psi$. Our approach to modelling the phase variations at different wavelengths provides a template for how to separate thermal emission from reflected light in spectrally resolved *James Webb* Space Telescope phase curve data.

**Key words.** techniques: photometric – planets and satellites: atmospheres – planets and satellites: physical evolution – planets and satellites: individual: MASCARA-1 b

## 1. Introduction

Stars with types earlier than ~F6 usually rotate rapidly, with O-type stars being observed to reach line of sight projected

rotational velocities ($v \sin i_*$) as high as 610 km s$^{-1}$ (Ramírez-Agudelo et al. 2013). Unlike the convective zones present in the outer layer of later types, the radiative outer layers of these stars lack the magnetic activity that dissipates their angular momentum with time, with the rapid rotation causing the radius at the equator to bulge with respect to the poles. von Zeipel (1924) predicted that an effect of this oblateness would be that the flux

<footnote>* The photometric time series data are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/658/A75</footnote>
emitted from the surface of the star would vary with local effec-
tive gravity. This results in both the brightness and temperature
decreasing from the poles to the equator: a phenomenon known
as gravity darkening.

The combined effect of the oblateness and surface tempera-
ture variations in the star causes the light curves of any transiting
planet to deviate from the characteristic transit signatures of
cooler hosts significantly (see Barnes 2009). Fitting for these
deviations can reveal a wealth of new information about the
system including the spin-orbit alignment, which provides an
insight into the formation and evolution of the system.

Spectroscopy of exoplanets during transit has proven to be a highly successful method of measuring their spin-orbit alignment. Observations of the Rossiter-McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924) have been the most productive source of these measurements to date for slow-rotating stars (e.g. Queloz et al. 2000; Brothwell et al. 2014; Hirano et al. 2020). Although these observations facilitate a measurement of the sky-projected spin-orbit angle $\lambda$, the fact that the stellar inclination $i_\star$ is degenerate with the stellar rotational velocity $v$ means that a precise measurement of the true, three-dimensional spin-orbit angle $\Psi$ is not possible for the majority of cases where $v$ is difficult to measure. Doppler tomography observations can yield a similar, reliable measurement for faster rotators with few spectral lines (e.g. Collier Cameron et al. 2010; Gaudi et al. 2017; Dorval et al. 2020), but studies to date have not been able to resolve $i_\star$. A common method of breaking the aforementioned degeneracy is the use of asteroseismology to estimate $i_\star$ (e.g. Chaplin et al. 2013; Benomar et al. 2014; Stefansson et al. 2020), although analytical criteria derived by Kamiaka et al. (2018) suggest that this method is only reliable for the range of $20^\circ < i_\star < 80^\circ$. Another method is the reloaded RM technique defined in Cegla et al. (2016), although the requirements of radial velocity measurements and a high-precision transit light curve has so far limited its application to a small number of systems (e.g. Bourrier et al. 2018; Dalal et al. 2019; Casasayas-Barris et al. 2021).

Unlike most spectroscopic methods, the transit light curve of a gravity-darkened host simultaneously encodes information about $v$, $i_\star$, and $\lambda$. This facilitates a measurement of $\Psi$ on the basis of time-series photometry alone, making it a powerful tool to study the history of planets orbiting early-type stars. Due to the high photometric precision required, studies have been limited until recently to planets and low-mass stellar companions orbiting a handful of the brightest stars in the *Kepler* field (Barnes et al. 2011, 2015; Szabó et al. 2011, 2020; Masuda 2015; Zhou & Huang 2013; Ahlers et al. 2014, 2015, 2019; Howarth & Morello 2017; Herman et al. 2018). Observations from the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) have the benefit of observing many more of the systems producing the most pronounced gravity darkening signatures than *Kepler*. However, its smaller collecting area, shorter observing baselines, and redder and narrower bandpass limit spin-orbit measurements to particularly favourable systems (e.g. Szabó et al. 2020; Ahlers et al. 2020b,a). To date, studies of this nature have only been conducted at optical wavelengths, and all but one (Szabó et al. 2020) have only used data taken at one wavelength band.

The main drawback of the gravity darkening method is the extent to which many of the observables from the transit light curve are correlated with each other. This makes accurate measurement of the full range of observables very challenging, and results are often highly sensitive to the choice of fixed parameters and the placement of priors within the model. The clearest example of this is the case of Kepler-13A b, where measurements

of $\lambda$ range from a slightly misaligned value of 23 ± 4 deg (Barnes et al. 2011) to a near-polar value of 59.20 ± 0.05 deg (Howarth & Morello 2017): the latter of which is in close agreement with a measurement from Doppler tomography (Johnson et al. 2014).

Roughly half an orbital period after the transit, a smaller drop in brightness is observed for the duration that the dayside of the exoplanet is occulted by its host. Optical wavelengths are an important window in which to observe the occultations of hot Jupiters, as the sharp drop off in thermal emission in the optical allows the reflective properties of their daysides to be studied. With the exception of Kepler-7b (e.g. Demory et al. 2013), the results from the *Kepler* and K2 missions suggested that the population of observed hot Jupiters reflect a tiny fraction of the starlight incident upon their daysides (systematic studies of *Kepler* phase curves are presented in Angerhausen et al. 2015; Esteves et al. 2015). Results from TESS satellite support this, though Wong et al. (2020a) and Wong et al. (2021) observe a tentative trend of increasing $A_g$ with planetary equilibrium temperature. Evans et al. (2013) use *Hubble/STIS* to report a significant detection of the occultation of HD 189733 b at blue wavelengths – in contrast to a non-detection at redder wavelengths – perhaps a signature of Rayleigh scattering. The limited results from ground-based occultation observations in the optical and near-ultraviolet include similar non-detections (e.g. Chen et al. 2014; Hooton et al. 2018), but also disagreeing results that are challenging to interpret (a full discussion is contained in Hooton et al. 2019). Precise occultation observations of a greater number of hot Jupiters at multiple wavelengths are required to further understand the atmospheric processes shaping the reflective properties of these extreme planets.

The Characterising Exoplanet Satellite (CHEOPS; Benz et al. 2021) was launched in December 2019, and has been acquiring high-precision optical photometry since going into full scientific operation in April 2020. Unlike the observing strategies of *Kepler* and TESS, CHEOPS is a mission designed primarily to observe individual targets already known to host transiting planets, with the first results being presented in literature in recent months (Bonfanti et al. 2021; Borsato et al. 2021; Delrez et al. 2021; Leleu et al. 2021; Morris et al. 2021; Swayne et al. 2021; Szabó et al. 2021; Van Grootel et al. 2021). With photometric precision and a wide optical bandpass<fnref n="1" /> both comparable to that of *Kepler*, this strategy allows CHEOPS to extend photometric measurements of $\Psi$ and occultation depth to the most amenable targets across most of the night sky. Demonstrating this capability, Lendl et al. (2020) reveal the asymmetry in the transit light curve of WASP-189 b at high precision, confirming the polar $\lambda$ reported by Anderson et al. (2018) and measuring $\Psi = 86.4_{-4.4}^{+2.9}$ deg. In addition, the occultation depth of 87.9 ± 4.3 ppm measured from four CHEOPS occultation light curves suggests a geometric albedo $A_g$ consistent with 0, although with no existing occultation data at longer wavelengths to disentangle thermal and reflected components of flux, an $A_g$ as high as ~0.3 is also possible.

The Jupiter-sized MASCARA-1 b (Talens et al. 2017, T17 hereafter) orbits its A8 type host star HD 201585 (MASCARA-1 hereafter) with a period of ~2.14 days. The resulting planetary equilibrium temperature in excess of 2500 K places it squarely amongst the hottest, most highly irradiated known exoplanets. Doppler tomography of a MASCARA-1 b transit also in T17 reveals both that the star is a rapid rotator

<footnote n="1">¹ The CHEOPS bandpass is almost identical to Gaia $G_{mag}$. More information can be accessed at https://www.cosmos.esa.int/web/cheops/performances-bandpass</footnote>
$(v \sin i_{\star} = 109 \pm 4 \text{ km s}^{-1})$, and that its planetary companion orbits in a near-polar orbit ($\lambda = 69.5 \pm 3.0^{\circ}$). With Gaia $G_{\text{mag}} = 8.25$, it is one of the brightest stars known to host a hot Jupiter, making it particularly amenable to classification using high-precision space-based photometry. These aspects are all analogous to the well-studied ultra-hot Jupiters KELT-9 b (Gaudi et al. 2017) and WASP-33 b (Collier Cameron et al. 2010), albeit with a slightly larger separation from its host in each case. Recently, Bell et al. (2021) present the analysis of a *Spitzer/IRAC* 4.5 µm full-phase curve of MASCARA-1 b acquired during its final year of operation as part of a population survey of *Spitzer* phase curves. However, MASCARA-1 eluded observation during the *Kepler* and K2 missions, and has not yet been observed by TESS due to its declination close to the ecliptic plane.

In this paper, we present Spi-OPS: the first joint analysis of *Spitzer/IRAC* and CHEOPS photometry. In it we present new CHEOPS observations of the asymmetric transit and occultation of MASCARA-1 b, which we analysed jointly with the *Spitzer* phase curve. Using this, we characterised the oblate star, the misaligned planetary orbit, and the atmosphere of the highly irradiated exoplanet. We also include an investigation into the effect of using transit photometry at different wavelengths in measuring properties such as stellar rotation, stellar inclination and spin-orbit alignment. We present an independent derivation of the stellar parameters and discuss the possibility that the star is pulsating in Sect. 2; we describe the observations and data reduction in Sect. 3; we describe how the light curves were jointly modelled Sect. 4; we investigate the effect of the placement of different priors when fitting CHEOPS and *Spitzer* transits in Sect. 5; we search for evidence of orbital precession in Sect. 6; we discuss the results of the joint model and present modelling of the atmosphere of the dayside in Sect. 7, and conclude in Sect. 8.

## 2. Stellar classification

### 2.1. Derivation of stellar parameters

In the following section, we describe our independent derivation of various stellar parameters including effective temperature $T_{\text{eff}}$, average radius $R_{\star}$, mass $M_{\star}$, and age $t_{\star}$.

We used a HARPS high resolution spectra with a signal-to-noise ratio (S/N) of ~1000 at the central wavelength of 550 nm to derive $T_{\text{eff}}$ from the broad line wings of H$\alpha$ with the publicly available spectral analysis package Spectroscopy Made Easy (SME; Valenti & Piskunov 1996; Piskunov & Valenti 2017). We used SME version 5.2.2, selected the Atlas12 (Kurucz 2013) stellar atmosphere grids, and retrieved the line data from VALD (Ryabchikova et al. 2015). We fixed the projected stellar rotation velocity to 109 km s$^{-1}$ from T17 and obtained $T_{\text{eff}} = 7490 \pm 150$ K, which is in agreement with McDonald et al. (2012) and T17. The iron abundance relative to hydrogen was modelled from iron lines in the range 6200–6800 Å, and was found to be $[\text{Fe/H}] = 0.15 \pm 0.15$.

We determined $R_{\star}$ using the infrared flux method (IRFM Blackwell & Shallis 1977) via relationships between the bolometric flux, the stellar angular diameter, the effective temperature, and the parallax. Using a Markov-chain Monte Carlo (MCMC) algorithm, we used the stellar spectral parameters derived above as priors to build spectral energy distributions from the ATLAS Catalogues (Castelli & Kurucz 2003). We subsequently computed synthetic fluxes by convolving these models with the throughput of the considered photometric bands that were compared with the observed fluxes in these bandpasses;

<table><caption>**Table 1.** Properties of MASCARA-1 and the methods employed to derive them.</caption><thead><tr><th colspan="3">MASCARA-1</th></tr></thead><tbody><tr><td>Alternative<br>names</td><td colspan="2">HD 201585<br>HIP 104513<br>Gaia DR2 1744911763437512064</td></tr><tr><th>Parameter</th><th>Value</th><th>Method</th></tr><tr><td>V (mag)</td><td>8.27</td><td>Simbad</td></tr><tr><td>G (mag)</td><td>8.2424</td><td>Gaia EDR3</td></tr><tr><td>J (mag)</td><td>7.819</td><td>2MASS</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>7490 ± 150</td><td>spectroscopy</td></tr><tr><td>[Fe/H] (dex)</td><td>0.15 ± 0.15</td><td>spectroscopy</td></tr><tr><td>R<sub>⋆</sub> (R<sub>☉</sub>)</td><td>2.072 ± 0.022</td><td>IRFM</td></tr><tr><td>M<sub>⋆</sub> (M<sub>☉</sub>)</td><td>1.825<sup>+0.097</sup><sub>-0.095</sub></td><td>isochrones</td></tr><tr><td>t<sub>⋆</sub> (Gyr)</td><td>0.8 ± 0.2</td><td>isochrones</td></tr><tr><td>ρ<sub>⋆</sub> (ρ<sub>☉</sub>)</td><td>0.205 ± 0.013</td><td>from R<sub>⋆</sub> and M<sub>⋆</sub></td></tr><tr><td>L<sub>⋆</sub> (L<sub>☉</sub>)</td><td>12.1 ± 1.0</td><td>from R<sub>⋆</sub> and T<sub>eff</sub></td></tr><tr><td>R<sub>⋆</sub> (R<sub>☉</sub>)</td><td>2.1 ± 0.2</td><td>T17</td></tr><tr><td>M<sub>⋆</sub> (M<sub>☉</sub>)</td><td>1.72 ± 0.07</td><td>T17</td></tr><tr><td>ρ<sub>⋆</sub> (ρ<sub>☉</sub>)</td><td>0.23<sup>+0.03</sup><sub>-0.01</sub></td><td>T17</td></tr></tbody></table>

Gaia G, $G_{\text{BP}}$, and $G_{\text{RP}}$ (Gaia Collaboration 2021); 2MASS J, H, and K (Skrutskie et al. 2006); and WISE W1 and W2 (Wright et al. 2010) to obtain the stellar bolometric flux and hence the angular diameter. Combined with the offset-corrected Gaia EDR3 parallax (Lindegren et al. 2021), we determined the stellar radius to be $2.072 \pm 0.022 R_{\odot}$.

We inferred $M_{\star}$ and $t_{\star}$ from two different sets of stellar evolutionary models, namely PARSEC v1.2S (Marigo et al. 2017) and CLÉS (Scuflaire et al. 2008), by adopting $T_{\text{eff}}$, [Fe/H], and $R_{\star}$ as input parameters. The output values from PARSEC v1.2S were computed interpolating within the grids of models by applying the Isochrone placement algorithm presented in Bonfanti et al. (2015, 2016). In the case of CLÉS, a direct fit to the stellar models was instead performed. After checking the consistency of the two pairs of parameters through the validation procedure described in detail in Bonfanti et al. (2021), we combined the probability distributions of both $M_{\star}$ and $t_{\star}$ to retrieve the corresponding medians and standard deviation as our reference values. A selection of parameters derived through the processes described in this section are shown in Table 1.

### 2.2. Possible pulsations

The fundamental parameters of MASCARA-1 mean that it is likely to exhibit one of or both $\delta$ Scuti- and $\gamma$ Dor-type stellar pulsations. The former class corresponds to pulsations with periods typically between half an hour to several hours, with the latter between one and three days. The brightness variations induced by these pulsations are mostly the consequence of variations in the effective temperature. Assuming a black-body radiation spectrum, the largest brightness-variation amplitudes would be detected at wavelengths of ~400 nm (covered at the blue end of the CHEOPS bandpass) for a star like MASCARA-1, while no detectable variation is expected in the infrared domain of *Spitzer*.

We checked which pulsation modes are predicted to be unstable for the MASCARA-1 stellar model computed with the CLÉS code. We used the non-adiabatic oscillation code MAD (Dupret 2001) and found that only radial and non-radial modes of $\delta$ Scuti
<table><caption>**Table 2.** Details of each of the five CHEOPS observations of MASCARA-1.</caption><thead><tr><th>Filekey <sup>(a)</sup></th><th>Start time<br>(UTC)</th><th>End time<br>(UTC)</th><th>N<sub>exp</sub>/N<sub>int</sub> <sup>(b)</sup></th><th>t<sub>int</sub><br>(s)</th><th>Type</th><th>N<sub>frames</sub></th><th>Efficiency<br>(%)</th></tr></thead><tbody><tr><td>PR100016_TG008501_V0200</td><td>2020-07-01 12:34</td><td>2020-07-01 22:20</td><td>1</td><td>23.0</td><td>Occultation</td><td>897</td><td>58.7</td></tr><tr><td>PR100020_TG000401_V0200</td><td>2020-07-09 00:24</td><td>2020-07-09 12:36</td><td>2</td><td>43.6</td><td>Transit</td><td>577</td><td>57.2</td></tr><tr><td>PR100016_TG008502_V0200</td><td>2020-07-14 10:48</td><td>2020-07-14 20:17</td><td>1</td><td>23.0</td><td>Occultation</td><td>919</td><td>61.9</td></tr><tr><td>PR100016_TG008503_V0200</td><td>2020-07-23 01:08</td><td>2020-07-23 10:54</td><td>1</td><td>23.0</td><td>Occultation</td><td>1201</td><td>78.5</td></tr><tr><td>PR100020_TG001101_V0200</td><td>2020-07-24 01:27</td><td>2020-07-24 13:58</td><td>2</td><td>43.6</td><td>Transit</td><td>860</td><td>83.3</td></tr></tbody></table>

Notes. <sup>(a)</sup>The filekey is the unique identifier associated with each dataset processed by the DRP. <sup>(b)</sup>N<sub>exp</sub>/N<sub>int</sub> denotes the image stacking order.

type are predicted to be excited by the κ mechanism (Pamyatnykh 1999) in this model: the predicted pulsation periods range from ~0.50 to 2.53 h.

We did not observe strong evidence for sizable δ Scuti-like pulsations in the CHEOPS light curves described in Sect. 3.1, although at the high end of the period range each light curve only covers a few periods of the expected signal. Moreover, the 98.77 minutes duration of an orbit of the CHEOPS satellite falls within this range, meaning that known instrumental signals occur on this timescale. Therefore, we made no attempt to fit for any low-amplitude δ Scuti pulsations in the CHEOPS light curve, and observations with a much longer baseline would be needed to reveal any pulsations emanating from MASCARA-1 b. We did observe time-correlated trends in some of the CHEOPS light curves that were not removed by conventional methods for detrending CHEOPS data, albeit over much longer timescales than the predicted periods of δ Scuti pulsations. We discuss how we modelled these in Sect. 4.3.

# 3. Observations and data reduction

## 3.1. CHEOPS

The CHEOPS satellite hosts an f/8 Ritchey-Chrétien on-axis telescope with an effective diameter of ~30 cm, with a single frame-transfer back-side illuminated charge-coupled device (CCD) detector. CHEOPS observed two transits and three occultations of MASCARA-1 b in July 2020 as part of the Guaranteed Time Observation (GTO) programme, which are summarised in Table 2. Due to the limited downlink bandwidth of the CHEOPS satellite, observations with exposure time $t_{exp} < 22.65$ s have their frames stacked prior to downlinking<fnref n="2" />. The transit and occultation observations were acquired as part of two different GTO sub-programmes. Although very similar and chosen to achieve ~90% of pixel full-well capacity, the difference in requested $t_{exp}$ led to them coincidentally falling either side of the 22.65 s threshold, resulting in the transit observations having a stacking order of two and the occultations a stacking order of one. The result of this is a large difference in integration time $t_{int}$ (the total $t_{exp}$ per downlinked frame), visible in Table 2. As CHEOPS occupies a low-Earth orbit with a period of ~98.77 min, the listed efficiencies refer to the proportion of the allocated time that the target was observable. Interruptions are primarily caused by Earth occultations, the satellite was crossing the South Atlantic Anomaly (SAA), and the detected stray light exceeding a predefined threshold computed by the instrument team, which is designed to keep the CHEOPS noise budget within requirements.

Following downlink, the data were reduced using version 13 of the CHEOPS Data Reduction Pipeline (DRP; described in full in Hoyer et al. 2020). In short, each frame is calibrated (by applying corrections for bias, gain, non-linearity, dark current and flat fielding), corrected for environmental effects (cosmic rays, smearing trails from nearby stars, and background), and then aperture photometry is performed for three fixed aperture sizes, along with a fourth that is performed at a size that minimises scatter due to contamination from background sources. The DRP estimates the contamination from nearby objects by simulating the CHEOPS field of view for each frame in the observation based on the Gaia DR2 star catalogue (Gaia Collaboration 2021) and a template of the extended CHEOPS point spread function (PSF). Due to the brightness of MASCARA-1, the mean contamination through all of the CHEOPS observations was estimated to be <0.05%. We selected the DEFAULT aperture of 25 pixels (25′′) in all cases, which minimised the route mean square (RMS) of each of the light curves, and discarded data points where the flux fell >5σ from the mean of the ten closest data points. The light curves are shown for the transits – visibly asymmetric due to the stellar rotation-induced gravity darkening – in Fig. 1 and for the second and third occultations in Fig. 2, with the models described in Sect. 4 also shown in each case. Collectively, the CHEOPS light curves cover 43.5% of full-phase coverage.

During one CHEOPS orbit of the first occultation observation, a problem with the guiding caused the target PSF to occupy a position on the detector roughly ten pixels away from its position for the other orbits (see top panel of Fig. 3). For bright targets like MASCARA-1, the CHEOPS payload is in-the-loop for improved pointing performances (jitter of order 1–2 arcsec). For the orbit in question, the centroiding algorithm could not lock the pointing on the target and the payload was consequently not in the loop. This issue occurred right after an interruption where the target was not visible for 40 minutes. Updates to the centroiding algorithm have improved its performance since then, and this issue has not been observed again.

When fitting an occultation model to this dataset alone, the measured occultation depth was very sensitive to the combination of auxiliary observing parameters used in the baseline model (see Sect. 4.3). This is because the erroneous orbit fell very close to the end of the observation, with the ten integrations in the final orbit proving insufficient to discriminate between a range of baseline models. This is demonstrated in the bottom panel of Fig. 3, which shows the raw data for the first occultation. The flux predicted by models that include detector position (orange) and time (violet) in the baseline model significantly diverge during the final few orbits, with the respective occultation depths measured to be 232 ± 19 and 164 ± 18 ppm. Not only do these significantly disagree with each other, but they are much

<footnote n="2">² See the CHEOPS observers manual.</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Normalised flux associated with the two CHEOPS transit observations as a function of orbital phase, with an arbitrary offset applied for display purposes. The raw data are shown in light blue points, and data in bins of 12 min are shown in dark blue points with error bars. The best fitting model is shown with a green line, and the 32 models evenly spaced in the chain are shown in fainter green lines. The *top panel* shows the raw flux, and the *bottom panel* shows the residual flux. The asymmetry visible in the transits is caused by stellar rotation-induced gravity darkening. The root mean square (RMS) of the residuals is 365 ppm.

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

**Fig. 2.** Similar to Fig. 1, but displaying the second and third CHEOPS occultation observations. The dashed lines mark the phases of the beginning and end of the occultation. The RMS of the residuals is 226 ppm.

deeper than the depths measured for the other two occultations, which are both relatively insensitive to the choice of linear basis vectors used. For this reason we excluded the first occultation from the joint analysis described in the following section. However, we cannot rule out the possibility of part of the increase in

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

**Fig. 3.** *Top: x* (indigo) and *y* (beige) positions of the PSF centroid on the detector. For the orbit that falls between the dotted lines, the PSF occupies a position roughly ten pixels away compared to the other orbits. *Bottom: similar to the top panel of Fig. 2*, but displaying the first CHEOPS occultation observation. The orange lines display models including *y* and *y*² detector position in the baseline model. The violet lines display models instead including time in the baseline model.

depth materialising due to an increase in the brightness from the MASCARA-1 b dayside.

### 3.2. Spitzer

A full phase curve of MASCARA-1 b (which was previously presented in Bell et al. 2021) was observed by *Spitzer/IRAC* using the channel 2 4.5 µm bandpass as part of programme 14059 with PI Jacob Bean, which we downloaded from the *Spitzer* Heritage Archive<fnref n="3" />. The observation began at 2019-03-03 22:32 UTC shortly before the ingress of a MASCARA-1 b occultation, and ended after 3.50 days of continuous observation and shortly after the egress of the following occultation. In this time, 106496 32x32 pixel images (with a 39'' × 39'' field of view) were acquired with an exposure time of 1.92 s, which were packaged in 1664 data cubes of 64 frames and acquired across three Astronomical Observing Requests (AORs). A full description of the method used to reduce the *Spitzer* data is given in Demory et al. (2016), which produces the raw light curve shown in the top panel of Fig. 4, with changes in AORs marked with dotted black lines.

## 4. Light curve fitting

A good model for jointly fitting the CHEOPS and *Spitzer* data has to take account of a number of signals present: the asymmetric transits of the oblate star, the flux deficit during the occultations, the phase-dependent modulation as the visible hemisphere of the tidally locked planet rotates with respect to the observer, the instrumental trends that are routinely present in CHEOPS photometry (e.g. Lendl et al. 2020; Bonfanti et al. 2021), and any other time-correlated trends. To incorporate each of these

<footnote n="3">³ http://sha.ipac.caltech.edu</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Normalised flux associated with the *Spitzer* phase curve observation. The raw data are shown in orange points, and data in bins of 30 min are shown in brown points with error bars. The best fitting model is shown with a red line, and the 32 models evenly spaced in the chain are shown in fainter red lines. The times when a new AOR begins are marked with dashed vertical lines. The top panel shows the raw flux, the middle panel shows the flux corrected for the instrumental effects with the transit, occultations and phase variations, and the bottom panel shows the residual flux. The mean absolute deviation of the residuals is 244 ppm.</figcaption></figure>

signals, we created a custom model to simultaneously fit the light curves based on the existing tools available in *pycheops*<fnref n="4" /> (see Maxted et al. 2021), a package that allows simultaneous modelling of transiting and occulting planets or eclipsing binaries along with the most common systematic trends present in CHEOPS datasets.

The model we used to fit the light curves is the sum of flux contributions from the star, planet, and instrumental systematics, where each term is normalised by the mean stellar flux when MASCARA-1 b is fully occulted by its host. We describe each of these terms in the following sections.

## 4.1. Stellar model

The flux from the star decreases during the transit of MASCARA-1 b, and is assumed to be constant out of transit. We modelled the gravity-darkened transits using the GravityDarkenedModel in *PyTransit*<fnref n="5" /> (Parviainen 2015), which implements the widely used model first presented in Barnes (2009)<fnref n="6" />. This will be described in full in Deline et al. (2021), but we describe the main steps here.

As described in Sect. 1, the surface temperature of a rapidly rotating star is well-described by the Maeder (2009) adaptation of von Zeipel’s theorem (von Zeipel 1924), which links the temperature and gravity at any point on the stellar surface by the equation. This is given by

$$ \frac{T(\theta)}{T_{\text{pole}}} = \left(\frac{g(\theta)}{g_{\text{pole}}}\right)^{\beta}, \qquad (1) $$

where $\theta$ is the stellar latitude, $g$ is local surface gravity, $g$ is polar surface gravity, and $\beta$ is the gravity darkening coefficient. *PyTransit* assumes that such a rapidly rotating star is well-approximated by a Roche model. The mass is concentrated at the centre of the star, the inner layers are spherical, and only the outermost layers of the star exhibit distortion due to the rotation. *PyTransit* computes the flux from the star in a grid of points on the surface projected onto the plane of the sky, and generates the light curve by evaluating the grid points that are hidden during a planetary transit as a function of time.

Changes in the light curve due to both the local temperature differences and oblateness induced by the rapid rotation are taken into account. Local temperature can be derived using Eq. (1), with $g$ determined by equating the gravitational and centrifugal forces:

$$ g(\theta) = -\frac{G M_{\star}}{r_{\theta}^{2}} \hat{u}_{r} + \left(\frac{2\pi}{P_{\text{rot}}}\right)^{2} R_{\star} \sin(\theta) \hat{u}_{\perp}. \qquad (2) $$

Here, $G$ is the gravitational constant, $M_{\star}$ is the stellar mass, $P_{\text{rot}}$ is the stellar rotational period, $R_{\star}$ is the stellar equatorial radius, and $\hat{u}_{r}$ and $\hat{u}_{\perp}$ are unit vectors radially from the centre and perpendicular to the stellar spin axis, respectively. $r_{\theta}$ is the latitude-dependent stellar radius accounting for the stellar oblateness $f_{\star}$, and is equal to $R_{\star}$ at the equator and $R_{\star}(1 - f_{\star})$ at the poles. In turn, $f_{\star}$ is given by

$$ f_{\star} = \frac{3\pi}{2G\rho_{\star}P_{\text{rot}}^{2}}, \qquad (3) $$

where $\rho_{\star}$ is the stellar density. The local flux is then derived from the local temperature by integrating the temperature-dependent stellar spectrum over the filter transmission function, taking into account the effects of limb darkening. As the effect of the gravity darkening varies as a function of wavelength, the light curve is

<footnote n="4">⁴ https://github.com/pmaxted/pycheops</footnote>

<footnote n="5">⁵ https://github.com/hpparvi/PyTransit</footnote>

<footnote n="6">⁶ We note an error in Barnes (2009), which is that all instances of $(1 - f^2)$ should be replaced with $(1 - f)^2$ in their Eq. 14. The current equation biases the coordinate of the oblate photosphere along the axis parallel to the line of sight.</footnote>
computed by integrating the stellar flux as a function of time over the filter response function. With $T(\theta)$ as an input PyTransit allows the use of a Planck function or a PHOENIX spectrum (Husser et al. 2013) to represent the stellar spectrum for each point on the surface. The PHOENIX spectrum is particularly desirable for the CHEOPS observations due to the large deviation in the spectra of early-type stars from that of a Planck function. As the spectra only extend to the reddest wavelengths in the CHEOPS bandpass, we used a Planck function for the Spitzer/IRAC 4.5 µm bandpass, which is much more representative of the star at infrared than at optical wavelengths. For the sake of computational efficiency, the stellar spectra were interpolated onto 15 wavelength bins for each response function. In addition to the model presented here, we obtained consistent results using a fit of only the CHEOPS transit light curves using TLCM<fnref n="7" /> (Csizmadia 2020) following the same methodology as Lendl et al. (2020).

The following terms in our parameterisation of the transit light curve are compatible with a symmetric transit, such as that of Mandel & Agol (2002). They are:

- the orbital period $P$,

- time of inferior conjunction $t_0$,

- transit depth $\delta_{\text{trans}} = (R_p/R_*)^2$,

- stellar density $\rho_*$,

- impact parameter $b$,

- eccentricity $e$ and argument of perirson $\omega_{\text{peri}}$, parameterised as $e \cos \omega_{\text{peri}}$ and $e \sin \omega_{\text{peri}}$,

- and the quadratic limb darkening parameters $q_1$ and $q_2$ (Kipping 2013).

A variety of different parameterisations of gravity-darkened light curves have been presented in literature, which has on occasion led to common names and symbols having different definitions between different papers. We use the same definitions and notation that are used by Masuda (2015), with the light curve asymmetry parameterised by the following additional terms.

- stellar equatorial radius $R_*$.

- sky-projected spin-orbit angle $\lambda$: the sky-projected angle between the orbital and stellar equatorial planes. This term is used interchangeably with (sky-projected) stellar obliquity by Masuda (2015) and by much of the related spectroscopy literature. We measure orbital plane anti-clockwise with respect to the stellar equatorial plane (sometimes termed $\alpha$), as opposed to the reverse direction (sometimes termed $\beta$).

- stellar inclination $i_*$: the angle between the stellar rotation axis and the observer's line of sight. This is related to the Barnes (2009) definition of stellar obliquity $\phi$ – the angle between the equatorial plane and the observer's line of sight – by the equation $i_* = 90^\circ - \phi$.

- stellar rotational velocity projected onto the observer's line of sight $v \sin i_*$,

- gravity darkening coefficient $\beta$,

- and stellar polar temperature $T_{\text{pole}}$.

This parameterisation differs from that of Masuda (2015) by the use of $v \sin i_*$ instead of the stellar rotation frequency $f_{\text{rot}}$, and $\rho_*$ instead of $M_*$: all of which are related by the equation

$$ \frac{v \sin i_*}{2\pi f_{\text{rot}}} = \left( \frac{3M_*}{4\pi\rho_*} \right)^{\frac{1}{3}} \sin i_* \quad (4) $$

Our choice of input parameters was informed by the desire to test the effect of placing Gaussian priors on observables from the Doppler tomography and our stellar characterisation (see Sect. 2).

As the Doppler tomography observation presented in T17 measures a prograde orbit at high significance, we placed a prior of $-90^\circ < \lambda < 90^\circ$ deg to exclude retrograde solutions. With this taken into account, the combination of the transit light curve and the Doppler tomography still cannot distinguish between solutions with $\{b, i_*, \lambda\}$ and $\{-b, -i_*, -\lambda\}$. By applying the constraint of $b > 0$ and letting $i_*$ vary in the full range of $-90^\circ < i_* < 90^\circ \text{ deg}^8$, our model will sample the only solution that exists within these limits. This is explained in greater detail in Sect. 5.2.

As described in Sect. 1, large degeneracies exist between the measurements of many of the parameters listed above, and previous studies (e.g. Barnes et al. 2011; Masuda 2015) have demonstrated that large disagreements in the measured asymmetric parameters can arise depending on the choice of free and fixed parameters within the model. We present a detailed investigation of this in Sect. 5, where we test the effect of fixing and varying different parameters when fitting the CHEOPS and Spitzer light curves separately and together.

## 4.2. Planetary model

The flux from the planet was modelled as the product of models to account for phase-dependent flux modulation and the occultation. We modelled the phase-dependent flux by self-consistently generating the components of flux due to thermal emission and reflected starlight. Although numerous phase curves of hot Jupiters have also exhibited signals associated with Doppler boosting and tidal ellipsoidal distortion (e.g. Shporer et al. 2014; Wong et al. 2020b; Owens et al. 2021), relations presented by Morris (1985), Morris & Naftilan (1993), and Shporer et al. (2019) suggest that the amplitude of these signals for MASCARA-1 would be roughly 10 ppm and 25 ppm, respectively. As we would be unable to significantly detect signals of these size in the Spitzer phase curve and due to the limited phase coverage of the CHEOPS data, we did not fit for these signals in our model.

To model the thermal phase variations in the CHEOPS and Spitzer light curves, we used a set of mathematical basis functions known as parabolic cylinder functions, which were derived by Heng & Workman (2014) to describe a heated fluid layer on a rotating sphere with frictional forces. Morris et al. (2022; M21 hereafter) demonstrate that these basis functions naturally describe the chevron-shaped feature in hot Jovian photospheres (Showman & Polvani 2011; Tsai et al. 2014) and are able to fit a sample of eight Spitzer phase curves with a small number of physically motivated parameters. A distinguishing feature of this approach is that it generates a single 2D temperature map, which is used to compute the phase-dependent flux in multiple bandpasses.

The full method is described in M21, but we give a brief summary of the steps here. Firstly, we generated a 2D temperature map as a function of planetary latitude $\theta$ and longitude $\phi$. This is given by the equation (1) of M21,

$$ T(\theta, \phi) = \bar{T} \left( 1 + \sum_{m,l}^{\ell_{\text{max}}} h_{m,l}(\theta, \phi) \right), \quad (5) $$

<footnote n="7">⁷ http://www.transits.hu/</footnote>

<footnote n="8">⁸ As for all of our runs without a Gaussian prior on $b$ we saw that $b \to 0$ in the sampling, but the Doppler tomography suggests that $b$ is nonzero with high significance, for the remainder of the study we assume that the solution with positive $b$ exists for $0 < i_* < 90^\circ \text{ deg}$. However, our results cannot distinguish between results with negative $b$ for any given solution.</footnote>
where $\bar{T}$ is the mean background temperature. The $h_{m,l}$ basis functions describe perturbations to the mean background temperature and are given by Eq. (258) of Heng & Workman (2014),

$$h_{m,l} = \frac{C_{m,l}}{\omega_{\text{drag}}^2 \alpha^4 + m^2} e^{-\bar{\mu}^2/2} [\mu m H_l \cos(m\phi) \\ + \alpha \omega_{\text{drag}} (2lH_{l-1} - \bar{\mu}H_l) \sin(m\phi)] \quad (6)$$

where $\alpha$ is a dimensionless fluid number that depends on the Rossby and Prandtl numbers and $\omega_{\text{drag}}$ is the dimensionless drag frequency. M21 demonstrated that synthetic temperature maps generated by general circulation models (GCMs) are well-fitted by $\alpha = 0.6$ and $\omega_{\text{drag}} = 4.5$; we hold these parameter values fixed as they control the latitudinal distribution of temperature, which is not constrained by real phase curves. Each $h_{m,l}$ mode is described by a pair of wavenumbers $(m, l)$, which are analogous to quantum numbers in the quantum harmonic oscillator; the power in each mode is quantified by $C_{m,l}$. Other quantities include $\mu = \cos\theta$ and $\bar{\mu} = \alpha\mu$. The physicist's Hermite polynomials are represented by $H_l$. M21 further showed that one needs to specify a hot spot offset $\Delta\phi$ within the fitting procedure that translates the entire temperature map in longitude.

After generating the temperature map, we then integrated the ratio of the planetary and stellar flux as a function of orbital phase $\xi$ (normalised such that transits occur at $\xi = \pm\pi$) according to Cowan & Agol (2011):

$$\frac{F_p}{F_*} = \frac{1}{\pi I_*} \left(\frac{R_p}{R_*}\right)^2 \int_0^{\pi} \int_{-\xi-\pi/2}^{-\xi+\pi/2} I_p(\theta, \phi) \cos(\phi + \xi) \sin^2(\theta) d\phi d\theta, \quad (7)$$

with planetary intensity $I_p$ given by

$$\int \lambda F_{\lambda} B_{\lambda}(T(\theta, \phi)) d\lambda \quad (8)$$

for a filter with photon flux response function $F_{\lambda}$ in units of electrons per photon, where $B_{\lambda}$ is a Planck function. For $I_{\star}$, we used a PHOENIX spectrum interpolated according to the values of $T_{\text{eff}}$, log $g$, and [Fe/H] listed in Table 1.

Alongside other parameters already defined in Sect. 4.1, for $l_{\text{max}} = 1^9$ the thermal phase modulation can be modelled using three parameters: $C_{1,1}^{10}$, $\Delta\phi$, and $\bar{T}$. Conceptually, these three parameters can be compared in an analogy to a sinusoidal phase curve model where $C_{1,1}$ is related to the semi-amplitude of the sinusoid, $\Delta\phi$ is a constant phase offset of the sinusoid, and $\bar{T}$ denotes the DC constant offset term from zero-flux. We allowed $\Delta\phi$ and $\bar{T}$ to vary uniformly, and $C_{1,1}$ to vary logarithmically. Given the temperature map yielded using these parameters, one may derive the Bond albedo $A_B$ and redistribution efficiency $\varepsilon$, which we describe in Sect. 7.3.

For the reflected light component of the phase variation, we used the general ab initio solutions for a semi-infinite atmosphere derived by Heng et al. (2021); we emphasise that ‘semi-infinite’ refers to the optical depth and not the spatial extent of the

atmosphere. Generally, the amplitude of a reflected-light phase curve is described by the geometric albedo $A_g$, while its shape is quantified by the integral phase function $\Psi_{\text{ph}}$. These quantities are in turn functions of the single-scattering albedo $\omega_{\text{scat}}$ and the scattering asymmetry factor $g$, which are fundamental scattering parameters. Forward, reverse, and isotropic scattering correspond to $g = 1, -1$ and 0, respectively. The widely used Henyey-Greenstein scattering phase function or reflection law (Henyey & Greenstein 1941) is employed. It is worth noting that this is a more general approach than assuming a Lambertian sphere, which is an artificial construct where an object is equally bright in reflected light regardless of viewing angle; see Dyudina et al. (2016) for multiple examples of how the reflected light phase curves of Solar System planets and moons are not well described by the Lambertian reflection law. Since no phase offset is ultimately detected, it is sufficient to use Eqs. (8) and (9) of Heng et al. (2021), which respectively describe $A_g$ and $\Psi_{\text{ph}}$ for a homogeneously reflective sphere. Upon obtaining $\Psi_{\text{ph}}$, it may be integrated over phase angle to obtain the phase integral $q$ (Russell 1916). The spherical albedo is then $A_s = qA_g$. Previously, the homogeneous-sphere solution was employed, in tandem with a double Henyey-Greenstein scattering phase function, to successfully fit the Cassini phase curves of Jupiter (Heng & Li 2021).

We used a single set of {$C_{1,1}$, $\bar{T}$, $\Delta\phi$, $\omega_{\text{scat}}$, $g$} to model phase variation in the light curves of both CHEOPS and Spitzer. The CHEOPS observations described in Sect. 3.1 cover well under half of an orbit of MASCARA-1 b across five separate visits, and thus only impose loose upper and lower limits on some parameters and do not meaningfully constrain others. The effect of fitting the phase curve model to the CHEOPS and Spitzer data simultaneously is to use the precise constraints availed by the full phase coverage of the Spitzer light curve to estimate the phase modulation in the CHEOPS light curves. Although some level of wavelength-dependence is expected for each of the terms, we assumed that the values for the atmospheric layer probed by the Spitzer observation are a good approximation for the layer probed by CHEOPS. This also requires us to neglect the role that chemical emission and absorption features may have on the thermal emission. Spectrally resolved observations in the future will be far better placed to detect such deviations than the light curves acquired using two broadband features that we are presenting.

We modelled the occultations using a Mandel and Agol transit model (Mandel & Agol 2002) with the limb darkening parameters fixed to 0. The model uses the same input parameters as the symmetric transit, along with $R_*$ to compute a delay in the occultation due to the difference in light travel time between the star and planet (roughly 20 seconds in the case of MASCARA-1). Occultation depth is not an input for the model, and is derived by recording the sum of the fluxes of the thermal and reflected components described above at the phase of the occultation centre.

### 4.3. Instrumental model

We modelled systematic trends in the CHEOPS light curves using a set of linear basis vectors of auxiliary observation parameters that are output by the DRP, such as sky background, contamination, detector smearing correction, detector xy-coordinates, and a ‘glint’ function (see Maxted et al., in prep.) created in pycheops by fitting a spline function to the out-of-transit or occultation flux as a function of telescope roll angle. We found that linear decorrelation only using a glint function per CHEOPS dataset and no other parameters optimised the

<footnote n="9">⁹ M21 showed that setting $l_{\text{max}} = 1$ provides a good approximation for the analysis of Spitzer phase curves. Here, only the Hermite polynomials $H_0 = 1$ and $H_1 = 2\bar{\mu}$ in Eq. (6) are used.</footnote>

<footnote n="10">¹⁰ For $l_{\text{max}} = 1$, $C_{1,1}$ is the only non-zero $C_{m,l}$ power coefficient, as $C_{m,0} \equiv 0$ and the perturbation arising from $C_{m,l}$ would be equal and opposite to that created by $-C_{m,-l}$ occupying the same value.</footnote>
Bayes factor (Kass & Raftery 1995, see the following section), with the exception of the final CHEOPS occultation where sky background was also included.

As discussed in Sect. 2.2, we observed low-frequency time-correlated trends in the CHEOPS light curves that we are unable to attribute to astronomical or instrumental sources. For this reason, we modelled them using a Gaussian process (GP) with a simple harmonic oscillator kernel using the celerite2¹¹ package (Foreman-Mackey et al. 2017; Foreman-Mackey 2018). Following its introduction by Gibson et al. (2012), GPs offer a flexible method to model correlated noise and accurately determine the parameters in transit light curves (Gibson 2014). The GP is parameterised by the frequency of the undamped oscillator $\omega_0$, $S_0$ is proportional to the power in the power spectral density at $\omega = \omega_0$, and the quality factor of the oscillator $Q$. We allowed the first two to vary on a logarithmic scale and fixed the latter to 10.

For each Spitzer AOR, we performed the same process using polynomial functions of x and y centroid position and x and y FWHM of the point response function up to the second order, including xy cross terms. The full set of linear basis vectors used for both the CHEOPS and Spitzer detrending is shown in Table A.1.

## 4.4. Posterior sampling

We fitted the trends present in the CHEOPS light curves using a pycheops Multivisit object. To derive the joint posterior distribution, we used the Dynamic Nested Sampler (see Skilling 2004, 2006; Higson et al. 2019) in the dynasty package (Speagle 2020). Unlike a Markov chain Monte Carlo (MCMC) algorithm (implemented in pycheops using the emcee package Foreman-Mackey et al. 2013), Nested Sampling estimates the Bayesian evidence $Z$ of the model, allowing for robust model comparison. We ran the fit with numerous combinations of the linear basis vectors described above, and selected the model that returned the highest estimate of $Z$.

Due to residual time-correlated noise in the Spitzer phase curve, for each of our fits we adopted the time-averaging method described in Winn et al. (2007). This involved determining the mean of the scaling parameter $\beta$ (not to be confused with the gravity darkening coefficient $\beta$) for bin sizes of 1–100 data points, and then running the sampling again with the uncertainties for each of the Spitzer data points boosted by a factor of $\beta$.

## 5. Accurate spin-orbit angle determination

Simultaneously acquiring accurate measurements of the parameters listed in Sect. 4.1 from a transit light curve presents a major challenge. Barnes (2009) described how the light curves spin-orbit aligned planets orbiting oblate stars mimic those of spherical stars at high impact parameters, which introduces a bias into the measurement of the transit depth. Additionally, correlations exist between so many of these parameters that an incorrect assumption about one can lead to biases in the measurements of several others.

The different reported values of sky-projected spin-orbit angle $\lambda$ for Kepler-13A b – where planet and host are comparable to MASCARA-1 b – provide a demonstration of this. Doppler tomography of the transit (Johnson et al. 2014) measures $\lambda = 58.6 \pm 2.0$ deg. For measurements using the Kepler photometry, this agrees well with Howarth & Morello (2017)

(59.20 ± 0.05 deg), but poorly with Barnes et al. (2011) (23 ± 4 deg) and Herman et al. (2018) (27.9 ± 1.1 deg<fnref n="12" />). A common feature of the latter two studies is that neither use priors from the Doppler tomography. In addition, they either fix the limb darkening coefficients or do not fit both components in the quadratic law. Masuda (2015) perform a detailed analysis using different combinations of priors and found that fitting both quadratic law limb darkening components could produce a result consistent with the tomography, and suggested using a Gaussian prior on $\lambda$ informed by the Doppler tomography.

Our goal in this section is to determine the best combination of uniform priors, Gaussian priors, and fixed parameters to allow accurate and informative measurements of certain parameters of interest. Of particular importance is the true spin-orbit angle $\Psi$, which as described in Sect. 1 is difficult to measure using other techniques. This is given by the relation:

$$ \cos \Psi = \cos i_* \cos i_p + \sin i_* \sin i_p \cos \lambda, \quad (9) $$

where $i_p$ is the orbital inclination. The asymmetric transit light curves of rapidly rotating stars contain information simultaneously about $i_*$, $i_p$, and $\lambda$. Thus, an accurate measurement of $\Psi$ will be possible in the event that these three parameters can be accurately determined.

### 5.1. Choice of priors

First of all, we considered what information was available for the MASCARA-1 system from which we could potentially construct relevant priors, whether these came from previous observations or simulations of this system and similar host stars. The stellar classification presented in Sect. 2 precisely estimates $R_*$, $M_*$ (and therefore $\rho_*$) and $T_{\text{eff}}$. We can infer the values of $q_1$ and $q_2$ from previous studies based on the stellar type and wavelength dependency (Claret & Bloemen 2011; Claret 2021). Similarly we can infer the value of $\beta$ (Claret 2016).

Finally, T17 have already observed the system during transit using Doppler tomography, which facilitated a precise measurement of the impact parameter $b$, sky-projected rotational velocity $v \sin i_*$, and $\lambda$. As these parameters are derived effectively by measuring the geometric path of the planet across its host during the observations, we consider these measurements to be particularly robust against systematic biases. Further, the comparison of our estimations of $b$, $v \sin i_*$, and $\lambda$ when they are left free to vary in the fit with those obtained from tomography is a useful indicator of accuracy. Placing a prior on $v \sin i_*$ is relatively common practice when fitting gravity-darkened light curves (e.g. Masuda 2015; Ahlers et al. 2020b). The same strategy has been adopted less commonly for $\lambda$ (e.g. Masuda 2015), but there are no examples of this strategy being adopted for $b$ (or equivalently $i_p$).

To investigate how the results extracted from fitting the asymmetric transit compared to those from the Doppler tomography, we performed a series of different fits to the CHEOPS and Spitzer transits separately and jointly. In our initial tests where we let all of the transit parameters vary uniformly, a particularly strong degeneracy existed between $i_*$ and $\beta$. This is because as $i_* \to 90$ deg, a smaller proportion of the hotter, polar regions of the star that give rise to the most pronounced transit asymmetry are visible, which mandates a larger value of $\beta$ to fit any asymmetry, and vice versa. Even when trying to mitigate against this by

<footnote n="12">¹² The published value is -27.9 ± 1.1 deg. We assume that the disagreement arises due to a definition of $\lambda$ in the opposite direction to the other studies.</footnote>

<footnote n="11">¹¹ https://github.com/exoplanet-dev/celerite2</footnote>
placing a Gaussian prior on $\beta$, the fits regularly returned results with $\beta$ at an unphysically high value. For this reason, we fixed $\beta = 0.199$ according to Claret (2016) for later runs.

We also placed a Gaussian prior $R_\star$ according the IRFM determination described in Sect. 2, as this value is poorly constrained by the model. In addition, we went one step further by fixing $T_{\text{pole}} = 7490$ K according to the spectral classification in Sect. 2, as it is also poorly constrained by the model, and is important in determining the phase variations as described in Sect. 4.2. We experimented with placing a prior on $\rho_\star$ according to our determinations of $R_\star$ and $M_\star$, but we let this vary uniformly as the light curves constrained it more strongly. Following the recommendation of Masuda (2015), we let the limb darkening coefficients $q_1$ and $q_2$ vary uniformly. Finally, we also let $\delta_{\text{trans}}$ and $i_\star$ vary uniformly within wide limits, and $P$ and $t_0$ within tight limits set by an initial fit to all of the transits. This choice of priors and fixed parameters are consistent for all of the fits presented in this paper, apart from the search for orbital precession performed in Sect. 6.

We then tested the effect of fitting different combinations of the CHEOPS and *Spitzer* transits, both with and without Gaussian priors informed by the tomography placed on *b*, $v \sin i_\star$, and $\lambda$. For each instrument combination we performed two fits: one with Gaussian priors on these three parameters, and another where they can vary uniformly within the limits described in Sect. 4.1. The results of fits for a selection of parameters are shown in Table 3, and corner plots displaying the correlations for the fit to the CHEOPS-only and CHEOPS + *Spitzer* light curves both with and without the Gaussian priors is shown in Figs. B.1 to B.4.

Some of the differences between the values for CHEOPS transits and *Spitzer* light curves with no priors from tomography are stark. $\delta_{\text{trans}}$ for CHEOPS is > 5% larger than that of *Spitzer*, which is an order of magnitude greater than is expected from wavelength-dependent atmospheric absorption for planets with large scale heights. This disagreement is likely driven by the more significant asymmetry in the CHEOPS light curves, which in turn favours a lower $i_\star$ where more of the hotter polar regions are visible to the observer. For the CHEOPS-only and *Spitzer*-only fits, the measured values are also highly sensitive to whether or not the priors are placed.

In contrast, our tests show that there are two approaches that result in good agreement between the sampled parameters. The first is the placement of the priors from tomography, where the results for the CHEOPS-only and CHEOPS+*Spitzer* fits are all within $2\sigma$ of each other. A larger disagreement is observed with the *Spitzer*-only fit, but this may be explained by the lower S/N for this transit. The second approach is to jointly fit the CHEOPS and the *Spitzer* transits together. In this case, the transits are forced to share the same transit parameters, which in turn greatly diminishes the degeneracies existing between them. While less precise, they are also in good agreement with those from the CHEOPS-only and CHEOPS+*Spitzer* fits where priors from tomography were used. Figures B.1 to B.4 display corner plots from both of the CHEOPS-only and CHEOPS+*Spitzer* fits, which demonstrates the effect of both of these approaches singularly and together. We discuss these two approaches in the following subsections.

## 5.2. Constraints from Doppler tomography

The joint consideration of time-series photometry and Doppler tomography has two major advantages. Firstly, they each rule out different orbital configurations that could otherwise be possible.

<table><caption>**Table 3.** Measurements of a selection of parameters for different combinations of transit light curves both with and without the placement of Gaussian priors on *b*, $v \sin i_\star$ and $\lambda$ informed by the Doppler tomography.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>no priors</th><th>Tomography priors</th></tr></thead><tbody><tr><td colspan="4"><strong>CHEOPS only</strong></td></tr><tr><td>δ<sub>trans</sub></td><td>%</td><td>0.6396<sup>+0.0043</sup><sub>-0.0053</sub></td><td>0.6192<sup>+0.004</sup><sub>-0.0048</sub></td></tr><tr><td>ρ<sub>★</sub></td><td>g cm<sup>-3</sup></td><td>0.3198<sup>+0.0025</sup><sub>-0.004</sub></td><td>0.2943<sup>+0.0028</sup><sub>-0.0031</sub></td></tr><tr><td>b</td><td></td><td>0.028<sup>+0.05</sup><sub>-0.022</sub></td><td>0.113 ± 0.012</td></tr><tr><td>q<sub>1</sub></td><td></td><td>0.153<sup>+0.034</sup><sub>-0.032</sub></td><td>0.227<sup>+0.026</sup><sub>-0.024</sub></td></tr><tr><td>q<sub>2</sub></td><td></td><td>0.126<sup>+0.088</sup><sub>-0.074</sub></td><td>0.45<sup>+0.054</sup><sub>-0.052</sub></td></tr><tr><td>i<sub>★</sub></td><td>deg</td><td>29.2<sup>+12.0</sup><sub>-8.2</sub></td><td>58.2 ± 2.7</td></tr><tr><td>λ</td><td>deg</td><td>-54 ± 26</td><td>-72.0<sup>+3.1</sup><sub>-2.9</sub></td></tr><tr><td>v sin i<sub>★</sub></td><td>km s<sup>-1</sup></td><td>68<sup>+37</sup><sub>-16</sub></td><td>105.7<sup>+3.6</sup><sub>-3.4</sub></td></tr><tr><td colspan="4"><strong>Spitzer only</strong></td></tr><tr><td>δ<sub>trans</sub></td><td>%</td><td>0.601<sup>+0.013</sup><sub>-0.021</sub></td><td>0.5928<sup>+0.0058</sup><sub>-0.006</sub></td></tr><tr><td>ρ<sub>★</sub></td><td>g cm<sup>-3</sup></td><td>0.303<sup>+0.01</sup><sub>-0.018</sub></td><td>0.2821<sup>+0.0033</sup><sub>-0.0034</sub></td></tr><tr><td>b</td><td></td><td>0.062<sup>+0.081</sup><sub>-0.046</sub></td><td>0.1192<sup>+0.0084</sup><sub>-0.0095</sub></td></tr><tr><td>q<sub>1</sub></td><td></td><td>0.0208<sup>+0.0051</sup><sub>-0.005</sub></td><td>0.0211<sup>+0.0048</sup><sub>-0.0046</sub></td></tr><tr><td>q<sub>2</sub></td><td></td><td>0.148<sup>+0.048</sup><sub>-0.049</sub></td><td>0.161 ± 0.05</td></tr><tr><td>i<sub>★</sub></td><td>deg</td><td>73<sup>+11</sup><sub>-20</sub></td><td>73.6<sup>+9.0</sup><sub>-9.3</sub></td></tr><tr><td>λ</td><td>deg</td><td>-36<sup>+52</sup><sub>-40</sub></td><td>-70.2 ± 3.0</td></tr><tr><td>v sin i<sub>★</sub></td><td>km s<sup>-1</sup></td><td>90<sup>+50</sup><sub>-61</sub></td><td>108.7<sup>+3.8</sup><sub>-3.9</sub></td></tr><tr><td colspan="4"><strong>CHEOPS+Spitzer</strong></td></tr><tr><td>δ<sub>trans</sub></td><td>%</td><td>0.605<sup>+0.011</sup><sub>-0.01</sub></td><td>0.6149<sup>+0.0035</sup><sub>-0.0043</sub></td></tr><tr><td>ρ<sub>★</sub></td><td>g cm<sup>-3</sup></td><td>0.2874<sup>+0.0086</sup><sub>-0.0087</sub></td><td>0.2889<sup>+0.0028</sup><sub>-0.0036</sub></td></tr><tr><td>b</td><td></td><td>0.017<sup>+0.03</sup><sub>-0.013</sub></td><td>0.110 ± 0.012</td></tr><tr><td>q<sub>1,C</sub></td><td></td><td>0.257<sup>+0.028</sup><sub>-0.025</sub></td><td>0.29<sup>+0.027</sup><sub>-0.026</sub></td></tr><tr><td>q<sub>2,C</sub></td><td></td><td>0.514 ± 0.086</td><td>0.383<sup>+0.044</sup><sub>-0.04</sub></td></tr><tr><td>q<sub>1,S</sub></td><td></td><td>0.0018<sup>+0.0034</sup><sub>-0.0013</sub></td><td>0.0119 ± 0.0053</td></tr><tr><td>q<sub>2,S</sub></td><td></td><td>0.46<sup>+0.35</sup><sub>-0.31</sub></td><td>0.148<sup>+0.05</sup><sub>-0.052</sub></td></tr><tr><td>i<sub>★</sub></td><td>deg</td><td>68.3<sup>+4.3</sup><sub>-6.1</sub></td><td>59.8<sup>+2.6</sup><sub>-2.3</sub></td></tr><tr><td>λ</td><td>deg</td><td>-69<sup>+13</sup><sub>-12</sub></td><td>-72.7<sup>+2.8</sup><sub>-3.0</sub></td></tr><tr><td>v sin i<sub>★</sub></td><td>km s<sup>-1</sup></td><td>125 ± 13</td><td>107.9<sup>+4.5</sup><sub>-3.3</sub></td></tr></tbody></table>

**Notes.** The correlations between parameters for the CHEOPS-only and CHEOPS+*Spitzer* fits are shown in Figs. B.1 to B.4.

This is demonstrated in the visualisation in the left panel of Fig. 6, which displays four different transit paths of MASCARA-1 b across its gravity-darkened host for the CHEOPS + *Spitzer* fit with tomographic priors. The shape of the light curve shows that MASCARA-1 b transits the hotter polar regions of the star in the first half of the transit (i.e. $-180 < \lambda < 0$ deg), which are described by the prograde solution marked in blue and the equivalent retrograde solution marked in yellow. The tomography confirms both that the orbit is prograde (i.e. $-90 < \lambda < 90$ deg), and that the majority of its path across the star is spent occulting red-shifted sections of the stellar disk (rotating away from the observer, the right half of the stellar disk), which are described by the solutions marked in blue and red. Only the blue solution collectively describes both the light curve and the tomography. For completeness, the right panel shows the equivalent transits
for the solutions with $-i_\star$ due to the symmetry of the star about the equatorial plane.

Secondly, in most circumstances Doppler tomography facilitates accurate determination of $b, v \sin i_\star$, and $\lambda$, but not $\Psi$ due to the lack of information about $i_\star$. Transit light curves of fast-rotating stars encode information about all of these parameters, but their accurate determination is more challenging. By placing priors informed by the tomography on these parameters, this allows a robust measurement of $i_\star$, which enables an accurate derivation of $\Psi$. We adopted this strategy for our joint fit to the full set of light curves in Sect. 7.

### 5.3. Constraints from multi-colour photometry

Our fits also suggest that in the absence of any Doppler tomography, accurate determination of transit parameters can also be achieved by the use of the CHEOPS and *Spitzer* data together. This allows us to report a measurement of $\Psi = 70.8_{-12.1}^{+10.8}$ deg independently from the Doppler tomography. Although the precision is much worse than when priors from tomography are used, the result is in close agreement with the measurement using priors and confidently rules out orbits without significant misalignment. It is possible that the use of multi-colour photometry in general – rather than the specific use of optical and mid-IR photometry – would also facilitate a similarly accurate result, but further study should be undertaken to confirm this.

To the best knowledge of the authors, this is the only study to date to attempt to use an IR transit light curve to measure the spin-orbit angle of a system. Due to the lack of clear gravity darkening-induced asymmetry in the *Spitzer* transit light curve, we tested the effect of decoupling the gravity darkening coefficient $\beta$ for each instrument and fixing it to 0 for *Spitzer*. This enabled us to test the performance of the model from Barnes (2009) compared to a symmetric model. The Bayes factor $\mathcal{B} = \frac{Z_{\text{Barnes}}}{Z_{\text{sym}}}$ for the two fits is 1.09, indicating that both fits perform similarly well but that the fit with the Barnes *Spitzer* transit is very marginally preferred. On this basis, we cannot claim to significantly detect gravity darkening-induced asymmetry in the *Spitzer* transit.

The biggest differences in the median models occur unsurprisingly during the transit, with a maximum difference of 250 ppm and an in-transit mean absolute deviation of 51 ppm: both well within the uncertainties of the *Spitzer* photometry. Additionally, with the exception of the *Spitzer* limb darkening coefficients, the fits returned values that were within $1\sigma$ of each other. As we have no compelling physically or data-driven reasons to suggest that Barnes model does not hold at IR wavelengths, we proceeded using a single value of $\beta$ to describe both the CHEOPS and *Spitzer* transits in the joint model presented in Sect. 7.

## 6. Orbital precession

We tested the possibility of using the transits to search for evidence of orbital precession, which would allow us to measure the stellar quadrupole moment $J_2$ of MASCARA-1 (Szabó et al. 2012). Although the only previous measurements using photometry (Szabó et al. 2012; Masuda 2015) leverage the four years of precise data acquired by *Kepler* (coupled with TESS data in the case of Szabó et al. 2020), the intervening 16 months between the CHEOPS and *Spitzer* transit observations presented an opportunity to see if significant precession had occurred in that time. In addition, we used the NITES $r'$-band transit from

<table><caption>**Table 4.** Measured impact parameter $b$ and transit duration $t_{14}$ for the individual fits to four transit light curves, as described in Sect. 6.</caption><thead><tr><th>Transit</th><th>b</th><th>t<sub>14</sub> (h)</th></tr></thead><tbody><tr><td>NITES</td><td>0.22<sup>+0.12</sup><sub>-0.14</sub></td><td>4.191<sup>+0.08</sup><sub>-0.13</sub></td></tr><tr><td><i>Spitzer</i></td><td>0.092<sup>+0.022</sup><sub>-0.028</sub></td><td>4.2698<sup>+0.008</sup><sub>-0.0084</sub></td></tr><tr><td>CHEOPS 1</td><td>0.1595<sup>+0.0087</sup><sub>-0.0089</sub></td><td>4.2381<sup>+0.0052</sup><sub>-0.0053</sub></td></tr><tr><td>CHEOPS 2</td><td>0.1036<sup>+0.0093</sup><sub>-0.011</sub></td><td>4.2654<sup>+0.0039</sup><sub>-0.0037</sub></td></tr></tbody></table>

2016 Sep 17 presented in T17¹³, increasing the time spanned by the observations to 46 months.

Precession manifests itself in transit light curves as a change in the impact parameter $b$ (alternatively $i_p$ or transit duration $t_{14}$) and $\lambda$. Previous measurements of $J_2$ for other hot Jupiters (e.g. Szabó et al. 2012; Masuda 2015; Watanabe et al. 2020; Borsa et al. 2021) are of the order $10^{-4}$. Adopting the same value for MASCARA-1, we would anticipate observing $db/dt \sim 0.002$ and $d\lambda/dt \sim 0.34$ deg yr$^{-1}$. The precision of the data will prevent us to detect variations at this low level, but $J_2$ is a poorly understood quantity and we cannot rule out the possibility that it is much larger for MASCARA-1.

In our first search, we fixed all transit parameters apart from $b$ and $\lambda$ to median values in the CHEOPS + *Spitzer* fit with prior informed by the Doppler tomography as described in the previous section, and let $b, \lambda$ and the baseline parameters be sampled uniformly. However, due to the lower significance of the asymmetry in the NITES and *Spitzer* transits, it was not possible to meaningfully constrain the corresponding measurement of $\lambda$ from these light curves. As all previous reports of observed precession have been detected on the basis of changing $b$ alone, our next search adopted the same strategy using a two-step process. Firstly, we fitted the transit light curves together. We fixed $b, v \sin i_\star, \lambda$ to the values measured by the Doppler tomography in T17, and let $P, t_0, \delta_{\text{trans}}$, and $\rho_\star$ and the limb darkening parameters per instrument ($q_1$ and $q_2$) vary uniformly. Secondly, we fitted each transit light curve individually, fixing the varying parameters to the median values in the previous step and allowing only $b$ and the baseline parameters to vary uniformly. The results are presented in Table 4.

The ground-based NITES transit favours a larger value of $b$ than the others, but was detected at relatively low significance and is within $1\sigma$ of the measurement of $b = 0.122 \pm 0.012$ (Talens et al. 2018) from the Doppler tomography acquired two weeks later. Although the measurements of $b$ from *Spitzer* and the second CHEOPS transit are in good agreement (and in turn in good agreement with the Doppler tomography), the first CHEOPS transit is $>2\sigma$ higher. This measurement would also appear unreliable, as we would not expect significant precession to occur in the two weeks that separated the two CHEOPS transit observations. As the ingress and egress were poorly sampled in the first CHEOPS transit due to data gaps in the observation, this resulted in the fit favouring a shorter transit duration with an associated larger $b$.

With the measurements from the NITES and first CHEOPS transits deemed unreliable, we found no evidence of the occurrence of significant orbital precession in the time between the

<footnote n="13">¹³ Although another NITES transit acquired on 2016 Jul 23 was also presented, we did not include this as the observation ended at the beginning of egress, making the detrending more susceptible to systematic bias.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Phase-folded *Spitzer* and CHEOPS photometry centred on the transits and occultations, using the same format as the previous light curves. Point with error bars show data in bins of 12 minutes. From top to bottom, the panels display the *Spitzer* transit, the CHEOPS transit, the *Spitzer* occultation, and the CHEOPS occultation, respectively.</figcaption></figure>

*Spitzer* and second CHEOPS transit. We discuss strategies for a possible future study to explore this in greater depth in Sect. 8.

# 7. Results and discussion

The CHEOPS transits, CHEOPS occultations, *Spitzer* phase curve, and detrended and phase-folded transits and occultations are displayed in Fig. 1, 2, 4, and 5, respectively, along with 32 models drawn from the posterior in each case. The results of the fit are displayed in Table 5 and A.1.

## 7.1. Spin-orbit alignment

As detailed in Sect. 5, in our transit only fits we obtained values of $v \sin i_\star$ and $\lambda$ consistent with Doppler tomography presented in T17 when they were allowed to vary freely. However, in the phase curve model we chose to place Gaussian priors on both parameters, due to the higher significance and the relative insensitivity to correlations with other parameters for the Doppler tomography measurement. Unsurprisingly, our measurements of $v \sin i_\star = 101.7_{-4.2}^{+3.5}$ km s$^{-1}$ and $\lambda = -69.2_{-3.4}^{+3.1}$ deg are both within $2\sigma$ of the Doppler tomography$^{14}$. Unlike the Doppler tomography, the transit light curve encodes information about the stellar inclination $i_\star$, which we measured to be $55.5_{-2.9}^{+2.3}$ deg. As this

breaks the degeneracy between $v$ and $i_\star$, we measured the stellar rotation period $P_{\text{rot}}$ to be $0.853_{-0.016}^{+0.017}$ days. This rapid rotation causes the stellar equator to bulge to a radius $104.39 \pm 0.18\%$ the size of the polar radius. By calculating the local surface gravity at both the equator and pole and utilising Von Zeipel's theorem (von Zeipel 1924), we measured a stellar equatorial effective temperature of $7357.8_{-4.9}^{+5.8}$ K, over 100 K lower than the temperature at the poles.

Using Eq. (9) we derived a true spin-orbit angle of $72.1_{-2.4}^{+2.5}$ deg. Our measurements of $\lambda$ (top panel) and $\Psi$ (bottom panel) are shown in Fig. 7 in the context of similar measurements collected in the TEPCat orbital obliquity catalogue$^{15}$ (Southworth 2011). Spectroscopic methods have facilitated many more measurements of $\lambda$ than $\Psi$, with stars where $T_{\text{eff}} < 7000$ K most likely to occupy aligned, prograde orbits. Stars with $T_{\text{eff}} > 7000$ K are less amenable to RM effect measurement so the sample is much smaller, but following the trend proposed by Winn et al. (2010), the results appear less clustered and more evenly distributed across the full range of prograde and retrograde orbits. MASCARA-1 b occupies a fairly sparsely populated section of the axes with a prograde near-polar orbit with misalignment and $T_{\text{eff}}$ comparable to measurements for Kepler-13Ab that agree with the Doppler tomography (Johnson et al. 2014; Masuda 2015; Howarth & Morello 2017). WASP-33 b, for example, appears at the edge of a larger cluster of planets with near-polar orbits of similar magnitudes, albeit with retrograde orbits.

<footnote n="14">14 Although the result of $\lambda = 69.5 \pm 3.0$ deg given by T17 is approximately the negative of our result, the difference arises due to their study using the $\beta$ definition of $\lambda$, while we adopted the opposite $\alpha$ definition (see Sect. 4.1).</footnote>

<footnote n="15">15 https://www.astro.keele.ac.uk/jkt/tepcat/obliquity.html</footnote>
<table><caption>**Table 5.** Priors and best-fitting values for the transit, occultation, and phase curve model parameters described in Sect. 4, along with a selection of derived parameters.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Prior</th><th>Best fit</th></tr></thead><tbody><tr><td colspan="4"><em>Transit, occultation and phase-curve model parameters</em></td></tr><tr><td>Period</td><td>P</td><td>d</td><td>U(2.1487, 2.1488)</td><td>2.14877381<sup>+0.0000088</sup><sub>-0.0000087</sub></td></tr><tr><td>Time of inferior conjunction</td><td>t<sub>0</sub></td><td>BJD<sub>TDB</sub> - 2458833</td><td>U(0.487, 0.489)</td><td>0.488151<sup>+0.000091</sup><sub>-0.000092</sub></td></tr><tr><td>Transit depth</td><td>δ<sub>trans</sub></td><td>%</td><td>U(0.5, 0.7)</td><td>0.6216<sup>+0.0035</sup><sub>-0.0033</sub></td></tr><tr><td>Stellar density</td><td>ρ<sub>★</sub></td><td>g cm<sup>-3</sup></td><td>U(0.2, 0.4)</td><td>0.2966<sup>+0.0027</sup><sub>-0.0024</sub></td></tr><tr><td>Impact parameter</td><td>b</td><td></td><td>U(0.0, 1.0)</td><td>0.113 ± 0.012</td></tr><tr><td>Cosine eccentric component</td><td>√e cos ω<sub>peri</sub></td><td></td><td>U(-0.15, 0.15)</td><td>0.001 ± 0.015</td></tr><tr><td>Sine eccentric component</td><td>√e sin ω<sub>peri</sub></td><td></td><td>U(-0.15, 0.15)</td><td>0.000<sup>+0.016</sup><sub>-0.015</sub></td></tr><tr><td>CHEOPS quadratic limb darkening</td><td>q<sub>1,C</sub></td><td></td><td>U(0.0, 1.0)</td><td>0.234<sup>+0.026</sup><sub>-0.022</sub></td></tr><tr><td>CHEOPS quadratic limb darkening</td><td>q<sub>2,C</sub></td><td></td><td>U(0.0, 1.0)</td><td>0.405<sup>+0.053</sup><sub>-0.049</sub></td></tr><tr><td>Spitzer quadratic limb darkening</td><td>q<sub>1,S</sub></td><td></td><td>U(0.0, 1.0)</td><td>0.00077<sup>+0.0018</sup><sub>-0.0006</sub></td></tr><tr><td>Spitzer quadratic limb darkening</td><td>q<sub>2,S</sub></td><td></td><td>U(0.0, 1.0)</td><td>0.46<sup>+0.35</sup><sub>-0.32</sub></td></tr><tr><td>Stellar radius</td><td>R<sub>★</sub></td><td>R<sub>☉</sub></td><td>N(2.072, 0.022)</td><td>2.082<sup>+0.022</sup><sub>-0.024</sub></td></tr><tr><td>Stellar inclination</td><td>i<sub>★</sub></td><td>deg</td><td>U(0, 90)</td><td>55.5<sup>+2.3</sup><sub>-2.9</sub></td></tr><tr><td>Sky-projected spin orbit angle</td><td>λ</td><td>deg</td><td>N(69.5, 3.0)</td><td>-69.2<sup>+3.1</sup><sub>-3.4</sub></td></tr><tr><td>Projected stellar rotational velocity</td><td>ν sin i<sub>★</sub></td><td>km s<sup>-1</sup></td><td>N(109.0, 4.0)</td><td>101.7<sup>+3.5</sup><sub>-4.2</sub></td></tr><tr><td>Gravity darkening coefficient</td><td>β</td><td></td><td>fixed</td><td>0.199</td></tr><tr><td>Stellar effective/polar temperature <sup>(a)</sup></td><td>T<sub>eff/pole</sub></td><td>K</td><td>fixed</td><td>7490</td></tr><tr><td>Hotspot offset</td><td>Δφ</td><td>deg</td><td>U(-180, 180)</td><td>2.0<sup>+8.9</sup><sub>-9.4</sub></td></tr><tr><td>C<sub>mf</sub> power coefficient</td><td>ln C<sub>1,1</sub></td><td></td><td>U(-2.75, 1.25)</td><td>-1.25<sup>+0.35</sup><sub>-0.38</sub></td></tr><tr><td>Mean background temperature</td><td>T̄</td><td>K</td><td>U(1000, 4000)</td><td>2350<sup>+130</sup><sub>-150</sub></td></tr><tr><td>Fluid number</td><td>α</td><td></td><td>fixed</td><td>0.6</td></tr><tr><td>Drag frequency</td><td>ω<sub>drag</sub></td><td></td><td>fixed</td><td>4.5</td></tr><tr><td>Single-scattering albedo</td><td>ω<sub>scat</sub></td><td></td><td>U(0.0, 1.0)</td><td>0.71<sup>+0.11</sup><sub>-0.18</sub></td></tr><tr><td>Scattering asymmetry factor</td><td>g</td><td></td><td>N(0.0, 0.07)</td><td>-0.010<sup>+0.070</sup><sub>-0.061</sub></td></tr><tr><td colspan="4"><em>Derived parameters</em></td></tr><tr><td>Stellar density</td><td>ρ<sub>★</sub></td><td>ρ<sub>☉</sub></td><td></td><td>0.2104<sup>+0.0019</sup><sub>-0.0017</sub></td></tr><tr><td>Stellar mass</td><td>M<sub>★</sub></td><td>M<sub>☉</sub></td><td></td><td>1.900<sup>+0.063</sup><sub>-0.068</sub></td></tr><tr><td>Transit duration</td><td>t<sub>14</sub></td><td>hr</td><td></td><td>4.226<sup>+0.010</sup><sub>-0.011</sub></td></tr><tr><td>Radius ratio</td><td>R<sub>p</sub>/R<sub>★</sub></td><td></td><td></td><td>0.07884<sup>+0.00022</sup><sub>-0.00021</sub></td></tr><tr><td>Semi major axis</td><td>a</td><td>R<sub>★</sub></td><td></td><td>4.1676<sup>+0.0047</sup><sub>-0.0051</sub></td></tr><tr><td>Semi major axis</td><td>a</td><td>au</td><td></td><td>0.040352<sup>+0.000046</sup><sub>-0.000049</sub></td></tr><tr><td>Orbital inclination</td><td>i<sub>p</sub></td><td>deg</td><td></td><td>88.45 ± 0.17</td></tr><tr><td>Eccentricity</td><td>e</td><td></td><td></td><td>0.00034<sup>+0.00034</sup><sub>-0.00023</sub></td></tr><tr><td>Argument of periastron</td><td>ω<sub>peri</sub></td><td>deg</td><td></td><td>-16<sup>+67</sup><sub>-52</sub></td></tr><tr><td>Planetary radius</td><td>R<sub>p</sub></td><td>R<sub>1</sub></td><td></td><td>1.597<sup>+0.018</sup><sub>-0.019</sub></td></tr><tr><td>Planetary equilibrium temperature</td><td>T<sub>equil</sub></td><td>K</td><td></td><td>2594.3<sup>+1.6</sup><sub>-1.5</sub></td></tr><tr><td>Stellar rotational period</td><td>P<sub>rot</sub></td><td>d</td><td></td><td>0.853<sup>+0.017</sup><sub>-0.016</sub></td></tr><tr><td>Oblateness</td><td>f<sub>obl</sub></td><td></td><td></td><td>0.0439 ± 0.0018</td></tr><tr><td>Stellar polar gravity</td><td>g<sub>pole</sub></td><td>m s<sup>-2</sup></td><td></td><td>125.6<sup>+1.7</sup><sub>-1.8</sub></td></tr><tr><td>Spin-orbit angle</td><td>Ψ</td><td>deg</td><td></td><td>72.1<sup>+2.5</sup><sub>-2.4</sub></td></tr><tr><td>CHEOPS occultation depth</td><td>δ<sub>occ,C</sub></td><td>ppm</td><td></td><td>133<sup>+19</sup><sub>-18</sub></td></tr><tr><td>Spitzer occultation depth</td><td>δ<sub>occ,S</sub></td><td>ppm</td><td></td><td>2116<sup>+53</sup><sub>-52</sub></td></tr><tr><td>Integrated dayside temperature</td><td>T<sub>d</sub></td><td>K</td><td></td><td>3062<sup>+66</sup><sub>-68</sub></td></tr><tr><td>CHEOPS nightside flux</td><td>F<sub>night,C</sub></td><td>ppm</td><td></td><td>0.96<sup>+2.60</sup><sub>-0.73</sub></td></tr><tr><td>Spitzer nightside flux</td><td>F<sub>night,S</sub></td><td>ppm</td><td></td><td>670<sup>+330</sup><sub>-310</sub></td></tr><tr><td>Integrated nightside temperature</td><td>T<sub>n</sub></td><td>K</td><td></td><td>1720 ± 330</td></tr><tr><td>CHEOPS phase curve amplitude</td><td>F<sub>amp,C</sub></td><td>ppm</td><td></td><td>131 ± 19</td></tr><tr><td>Spitzer phase curve amplitude</td><td>F<sub>amp,S</sub></td><td>ppm</td><td></td><td>1470<sup>+330</sup><sub>-340</sub></td></tr><tr><td>Geometric albedo</td><td>A<sub>g</sub></td><td></td><td></td><td>0.171<sup>+0.066</sup><sub>-0.068</sub></td></tr><tr><td>Phase integral</td><td>q</td><td></td><td></td><td>1.559<sup>+0.095</sup><sub>-0.100</sub></td></tr><tr><td>Spherical albedo</td><td>A<sub>s</sub></td><td></td><td></td><td>0.266<sup>+0.097</sup><sub>-0.100</sub></td></tr></tbody></table>

**Notes.** The baseline and noise parameters in the model are displayed in Table A.1. <sup>(a)</sup>One parameter is used for both the T<sub>eff</sub> (used in the phase-curve model) and T<sub>pole</sub> (used in the asymmetric transit model).
<figure><img src="image_7.png" /></figure>

**Fig. 6.** Visualisation of the path of the MASCARA-1 b transit across the stellar disk using the CHEOPS+Spitzer fit with priors from Doppler tomography. The star rotates anti-clockwise when viewed from the top, and the majority of the hemisphere in view moves from left to right. Constant stellar latitudes are marked with curved dashed lines. The colour scale is proportional to the flux emitted by the star, with brighter areas marked with lighter colours. For a given $i_*$, all possible transit paths of the planet across the star can be described using different combinations of $b$ and $\lambda$. The {$b, \lambda$} pair associated with the median path transited by the planet is shown with a solid black line with a blue arrow indicating the direction, and with 32 samples drawn from the posterior shown in grey. The dashed lines with red and orange arrows indicate alternative {$b, \lambda$} pairs that also describe the signals in the Doppler tomography and transit light curves, respectively. However, when considering the tomography and the transit light curve together, only the {$b, \lambda$} pair described by the blue arrow can successfully describe them both. Assuming stellar symmetry about the equatorial plane, there are two indistinguishable solutions for each arrow in the range $-90 < i_* < 90$ deg: one with {$b, i_*, \lambda$} and another with {$-b, -i_*, -\lambda$}. The {$b, \lambda$} pairs with positive $i_*$ (north pole pointing towards the observer) are shown in the left panel, and the {$b, \lambda$} pairs with negative $i_*$ (north pole pointing away from the observer) are shown in the right panel. The blue arrow in the left panel is an example of a solution fitting the light curve with $b > 0$, and $-90 < \lambda < 90$ deg, and corresponding solutions are found for all of the transit fits presented in this paper.

Recently, Albrecht et al. (2021) note that the majority of planets with misaligned orbits occupy nearly polar orbits in the range $\Psi = 80–120$ deg. Like Kepler-13A b, while relatively close to polar the $\Psi$ of MASCARA-1 b falls outside this particularly clustered range. Although not as massive as Kepler-13A b, with $M_p = 3.7 \pm 0.9 M_J$ measured by T17, MASCARA-1 b has a relatively large mass compared to other hot Jupiters, which may support a tentative trend of more massive misaligned planets being less likely to occupy polar orbits. Measurements for a much larger sample of planets transiting stars with $T_{\text{eff}} > 7000$ K will be necessary to reveal more detailed trends.

## 7.2. Phase curve analysis

As described in Sect. 4.2, the thermal component of the phase curve model is evaluated by generating a 2D temperature map, and integrating the blackbody flux between the respective filter response functions. The temperature map associated with our median result is displayed in Fig. 8. The chevron-shaped patterns of common local temperatures are the trademark signature of the $h_{m,l}$ basis functions, which replicate the shape of the global temperature distributions produced by many GCMs (e.g. Showman & Polvani 2011). As $\omega_{\text{scat}}$ is fixed to 4.5 in the fit, the hotspot offset $\Delta\phi$ is roughly equivalent to the offset between the phase of superior conjunction and maximum flux, which with $\Delta\phi = 2.0_{-9.4}^{+8.9}$ deg is consistent with zero.

A corner plot showing the correlations between the phase curve parameters is shown in Fig. 9. The negative correlation

between $\ln c_{1,1}$ and $\bar{T}$ is analogous to the relationship between the phase curve amplitude and the minimum flux from sinusoidal phase curve parameterisations. While the thermal flux variations are wavelength-dependent, the reflected light component contributes the same phase-dependent flux in all bandpasses. The nonzero values favoured in the sampling of $\omega_{\text{scat}}$, which scales the amplitude of the reflected flux, suggests that the CHEOPS occultation is too deep with respect to the *Spitzer* occultation to be described by blackbody flux alone.

Table 5 also lists a range of parameters derived from the phase curves. Although we list the CHEOPS phase curve parameters, we caution against interpreting any significance from the CHEOPS nightside flux value of $0.96_{-0.73}^{+2.60}$ ppm. This is derived from the posterior distributions of phase curve parameters fitting both the *Spitzer* and CHEOPS light curves simultaneously, but is predominantly informed by the *Spitzer* data.

Heng et al. (2021) describe how the geometric albedo $A_g$, spherical albedo $A_s$, and phase integral $q$ are derived from our parameterisation of the reflected light phase curve (see Sect. 4.2). Our measurement of $A_g = 0.171_{-0.068}^{+0.066}$ constitutes a $2.5\sigma$ detection of reflected light from the dayside (see Sect. 7.4 for a wider discussion). Our measurement of $q = 1.559_{-0.100}^{+0.095}$ is consistent with isotropic scattering ($q = 3/2$). Although the reflected light is only present in significant quantities in the CHEOPS bandpass and full-phase coverage is preferable when measuring $q$, the 43.5% full-phase coverage that the CHEOPS observations achieve are sufficient to meaningfully constrain this value.
<figure><img src="image_8.png" /></figure>

**Fig. 7.** All published measurements of the absolute value of the sky projected spin-orbit angle $|\lambda|$ (top panel) and true spin-orbit angle $\Psi$ (bottom panel) as a function of $T_{\text{eff}}$, according to TEPCat. All planets with measurements of $\Psi$ with $T_{\text{eff}} > 7000$ K are marked specifically. In cases where more than one measurement exists for a given planet, they are all shown. Some values of $\lambda$ were converted to fall in the range $-180 < \lambda < 180$ deg.

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

**Fig. 8.** Temperature map as a function of planetary latitude and longitude derived from the best fitting phase curve parameters. The colour bar on the right shows the corresponding temperatures. The longitude zero-point is defined with respect to the sub-stellar point.

These derived parameters facilitate a more direct comparison with the fit to the *Spitzer* phase curve previously presented in Bell et al. (2021). Our measurement of occultation depth and nightside flux – with the corresponding dayside and nightside integrated temperatures – are each $1-2\sigma$ higher than the Bell et al. (2021) equivalents. Although different approaches were

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

**Fig. 9.** Corner plot demonstrating the correlations between the five free parameters describing the phase variations. Increasing densities of data points are marked with contour lines and darker background shades. The top row displays a histogram showing the posterior distribution for the parameter associated with that column. Generated using the corner package (Foreman-Mackey 2016).

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

**Fig. 10.** Top panel: phase curve model associated with the median result for the *Spitzer* bandpass (maroon), which is the sum of the thermal and reflected light components described in Sect. 4.2. A Lambertian phase curve model with the same amplitude, minimum flux, and hotspot offset (green) is also shown. Bottom panel: flux difference between the two models.

employed between the two studies for light curve detrending, another likely provenance is the different phase curve parameterisations used, which is demonstrated in Fig. 10.

### 7.3. Energy budget: bond albedo and recirculation efficiency

Given the hemispherically averaged dayside and nightside temperature, one may derive the Bond albedo $A_B$ and redistribution factor $f$ using zero-dimensional (0D) ‘box models’ as stated in Eqs. (4) and (5) of Cowan & Agol (2011). However, this depends
on the choice of boundary condition associated with the night-
side temperature, for example demanding that *f* = 2/3 in the
limit of no heat redistribution. This issue may be avoided by
deriving the Bond albedo directly from the 2D temperature map,
as demonstrated in *Keating et al.* (2019) and M21.

$$
A_B = 1 - \left(\frac{a}{R_*}\right)^2 \frac{\int_{-\pi}^{\pi} \int_{0}^{\pi} F_p \sin \theta \, d\theta \, d\phi}{\pi \sigma_{SB} T_{\text{eff}}^4}, \quad (10)
$$

where σ<sub>SB</sub> is the Stefan-Boltzmann constant. The redistribu-
tion factor ε is simply the ratio of the flux from the integrated
nightside and dayside (Perna et al. 2012),

$$
\epsilon = \frac{\int_{\pi/2}^{\pi} \int_{0}^{\pi} F_p \sin \theta \, d\theta \, d\phi + \int_{-\pi}^{-\pi/2} \int_{0}^{\pi} F_p \sin \theta \, d\theta \, d\phi}{\int_{-\pi/2}^{\pi/2} \int_{0}^{\pi} F_p \sin \theta \, d\theta \, d\phi}. \quad (11)
$$

In the 2D approach, the redistribution factor *f* is undefined.

The derived values of ε and A<sub>B</sub> using both the 0D approach
from Cowan & Agol (2011) and 2D approach from M21 are dis-
played in Table 6, and in both cases they are within 1σ of each
other. Both methods yield a large proportion of samples with
negative values of A<sub>B</sub>, which suggests that it could be under-
estimated due to spectral features in the emission spectrum of
MASCARA-1 b within the *Spitzer*/IRAC 4.5 µm bandpass caus-
ing a significant deviation from the flux predicted by a Planck
function.

Formally, the Bond albedo is the spherical albedo weighted
by the stellar spectral energy distribution F<sub>⋆</sub>, integrated over all
wavelengths (Marley et al. 1999),

$$
A_B = \frac{\int_0^\infty A_s(\lambda) F_\star \, d\lambda}{\int_0^\infty F_\star \, d\lambda}. \quad (12)
$$

In this study, each quantity is derived separately: A<sub>B</sub> as described
in the previous paragraphs using the thermal component of the
light curve model and A<sub>s</sub> as the product of A<sub>g</sub> and the phase inte-
gral *q* from the reflected component of the light curve model.
Formally, the spherical albedo is defined at a single wavelength,
as is the case for the geometric albedo; both quantities are intrin-
sic properties of the atmosphere that do not formally involve
the star. In practice, spherical and geometric albedos are com-
monly reported as bandpass-integrated quantities, even for the
planets of the Solar System (e.g. Table 7 of Pearl & Conrath
1991); an exception is for Cassini data of Jupiter (e.g. Fig. 3 of Li
et al. 2018). Our CHEOPS bandpass-integrated<sup>16</sup> measurement
of A<sub>s</sub> = 0.266<sup>+0.097</sup><sub>-0.100</sub> is larger than our measurement of A<sub>B</sub> from
the 2D temperature map.

A possible cause of the relatively poor agreement between
these two values could arise due to the CHEOPS and
Spitzer/IRAC 4.5 µm bandpasses only covering broad opti-
cal and part of the mid-infrared, respectively. Consulting the
PHOENIX spectrum that modelled the stellar flux in the light
curve fit (shown in red in the bottom panels of Figs. 11 and
12), these bandpasses only capture ~69% of this flux. In the
event that MASCARA-1 b reflected none of the light incident
upon it at wavelengths outside of the CHEOPS and 4.5 µm band-
passes (i.e. A<sub>s,out</sub> = 0), and assuming that optical wavelengths

<table><caption>**Table 6.** Heat redistribution efficiency $\epsilon$ and Bond albedo $A_B$ derived using both the 0D method from Cowan & Agol (2011) and the 2D method using temperature maps from M21.</caption>
  <thead>
    <tr>
      <th>Parameter</th>
      <th>0D</th>
      <th>2D</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>&epsilon;</td>
      <td>0.23<sup>+0.13</sup><sub>-0.10</sub></td>
      <td>0.204<sup>+0.089</sup><sub>-0.068</sub></td>
    </tr>
    <tr>
      <td>A<sub>B</sub></td>
      <td>0.139<sup>+0.074</sup><sub>-0.100</sub></td>
      <td>0.057<sup>+0.083</sup><sub>-0.101</sub></td>
    </tr>
  </tbody>
</table>

and longer are well-described by the A<sub>s</sub> derived in our light
curve model (i.e. A<sub>s,in</sub> = 0.266<sup>+0.097</sup><sub>-0.100</sub>), using Eq. (12) we com-
puted A<sub>B</sub> = 0.184<sup>+0.067</sup><sub>-0.069</sub>. Although still over double the A<sub>B</sub> that
we derived for the temperature maps, they are only discrepant
at 1.2σ. The stellar flux that falls outside of these bandpasses is
split relatively evenly between the ultraviolet and near-infrared.
If MASCARA-1 b reflects a much lower proportion of incident
light from its host in the ultraviolet than in the optical, this would
be a reversal of reflected properties measured on the basis of
a Hubble/STIS secondary eclipse observation of HD 189733 b
(Evans et al. 2013). Reflection from hot Jupiters in the near-
infrared has not been well-studied to date, as the observed flux
in this range is generally dominated by thermal emission.

When assuming that the thermal emission from the dayside
follows a blackbody distribution (as our light curve model does),
the dayside integrated temperature of T<sub>d</sub> = 3062<sup>+66</sup><sub>-68</sub> K confirms
the position of the MASCARA-1 dayside amongst the hottest of
any known exoplanet. We report a 3σ detection of flux from the
nightside and a corresponding nightside integrated temperature
of T<sub>n</sub> = 1720 ± 330 K. When taken together with the hotspot
offset Δφ consistent with zero and the derived ε, this indicates
the relatively poor redistribution of heat from dayside to night-
side that is typical for the majority of hot Jupiters. However, a
comparison with the hot Jupiters analysed in Wong et al. (2021)
places the redistribution efficiency of MASCARA-1 b for both
the 0D and 2D methods above that of the majority of the sample.
Correspondingly, the extreme levels of irradiation MASCARA-
1 b receives from its host place the derived T<sub>n</sub> as the second
hottest in the sample, with only HAT-P-7 b hotter, although likely
also cooler than KELT-9b (Mansfield et al. 2020; Bell et al. 2021)
which did not form part of the sample. We add the caveat that
different analyses of Spitzer/IRAC phase curves have yielded
dramatically different conclusions about the energy budgets of
the corresponding exoplanets (see e.g. May & Stevenson 2020,
and references therein).

7.4. *Emission spectrum retrieval*

Although the previous section presented a self-consistent analy-
sis of the thermal and reflective properties, for the sake of com-
putational efficiency it necessarily included some major simpli-
fications. In particular, we assumed that the thermal emission
from each point in the 2D temperature map could be described
by a Planck function, and reduced the spectra for planet, star, and
filter transmission functions down to only 15 bins per bandpass.
To perform a more detailed spectral analysis of the dayside atmo-
sphere, we generated a grid of emission spectra using HELIOS<fnref n="17" />
(Malik et al. 2017, 2019) to jointly interpret the CHEOPS and
Spitzer occultation depths.

We used a planetary radius of 1.597 R<sub>j</sub> and an orbital dis-
tance of 0.040352 au derived from our light curve fitting, and

<footnote n="16">¹⁶ Although practically speaking we have actually measured A<sub>s</sub> collec-
tively across the CHEOPS and *Spitzer*/IRAC 4.5 µm bandpasses, in
reality the fact that the stellar flux in the 4.5 µm bandpass constitutes
~ 0.3% of the combined total means that it effectively approximates the
A<sub>s</sub> for the CHEOPS integrated bandpass only.</footnote>

<footnote n="17">¹⁷ https://github.com/exoclime/HELIOS</footnote>
<figure><img src="image_12.png" /></figure>

**Fig. 11.** *Top:* emission spectra generated to fit the *Spitzer* occultation depth, as described in Sect. 7.4. The spectrum from the HELIOS grid is shown in magenta and the best fitting blackbody is shown in orange. The unbinned spectra are shown with thin, lighter lines. Binned spectra are also shown in darker, thicker lines for display purposes. The measured CHEOPS and *Spitzer* occultation depths are shown with black squares with error bars, and the implied depths according to the spectra are shown with coloured circles (obscured by the measured depth in the case of *Spitzer*). The CHEOPS and *Spitzer* response functions are shown in grey. *Inset top:* the temperature–pressure profile of the best fitting HELIOS spectrum. *Bottom:* the same spectra, but shown in terms of flux instead of occultation depth. The PHOENIX spectrum used for the star is shown in red, shrunk by a factor of 50 for display purposes. Right inset: a closer look at the flux at wavelengths covered by the *Spitzer*/IRAC 4.5 µm bandpass. Other dotted lines in both panels are regions that are displayed in Fig. 12.

a surface gravity of $\log g = 3.55$ in cgs units informed by the planetary mass reported in T17. Since the metallicity of the host star is compatible with solar element abundances, for simplicity we correspondingly adopted these to describe the chemical composition of the planet. The considered opacity sources and corresponding references can be found in Wong et al. (2020b).

We generated a grid of self-consistent HELIOS atmosphere models for MASCARA-1 b as a function of the heat redistribution efficiency $\epsilon$ (see Lendl et al. 2020 for an example of such a grid constructed for WASP-189b). In a post-processing step, we then generated a high-resolution planetary emission spectrum without contributions by scattering for each point in the grid. These spectra, thus, contain only the thermal emission part of the planet's spectrum. To calculate the theoretical occultation depths implied by these emission spectra, we used the same PHOENIX spectrum as was used for the computation of the thermal flux modulation (see Sect. 4.2).

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

**Fig. 12.** Closer look at the emission spectra at wavelengths covered by the CHEOPS bandpass, enclosed by dotted lines in Fig. 11.

The high-resolution spectra are averaged over the *Spitzer* and CHEOPS bandpasses to obtain the corresponding, theoretical occultation depths as seen by these instruments. Since the bandpass-averaged occultation depths vary smoothly with the heat redistribution factor $\epsilon$, we parameterised them with a spline function for each filter.

We then used an MC algorithm to find the distribution of $\epsilon$ that best fit the *Spitzer* occultation depth within its error bars, and also taking into account the uncertainty on the transit depth and the $R_p/R_\star$ value reported in Table 5. Although flux originating from both thermal emission and reflected light comprised the model fitting the *Spitzer* occultation, this step implicitly assumes that all of the flux from the dayside in *Spitzer* bandpass arises due to thermal emission. As thermal emission is expected to dominate reflected light in the mid-infrared for ultra-hot Jupiters such as MASCARA-1 b (which is indeed true in our light curve model), we deem this to be an acceptable assumption.

This analysis yields a result of $\epsilon = 0.307 \pm 0.088$ ($f = 0.538 \pm 0.038$). The spectra for the median value of $\epsilon$ from this process are displayed in magenta in Fig. 11. As can be seen in the zoom of the spectra for the CHEOPS bandpass shown in Fig. 12, the CHEOPS occultation depth predicted using this method is $88 \pm 10$ ppm–significantly less than the observed value of $133_{-18}^{+19}$ ppm from our light curve model. Since we neglected scattering in the post-process emission spectra, the difference between the predicted occultation depth and the measured one can be translated into the geometric albedo $A_g$ in the CHEOPS bandpass, for which we obtained a value of $A_g = 0.133 \pm 0.062$.
We also tested the effect of replacing the emission spectrum described above with a blackbody (shown in orange in Figs. 11, 12), and fitting to the *Spitzer* occultation depth as a function of the brightness temperature $T_B$. This was optimised with $T_B = 3178 \pm 53$ K, which resulted in a predicted CHEOPS occultation depth of 90 ppm and an $A_g$ of $0.127 \pm 0.063$. The blackbody spectrum is broadly very similar to the HELIOS spectrum, with the main differences arising due to metal emission lines close to the *Spitzer* bandpass.

Given that the dayside emission spectrum produced by HELIOS is so closely approximated by a blackbody, the retrieval should return results in good agreement with those resulting from the light curve model, as is reported in the previous section. From this, we indeed see that the measurement of $A_g = 0.171_{-0.068}^{+0.066}$ is in good agreement with both the retrieved HELIOS and blackbody values, and $T_d = 3062_{-68}^{+66}$ K ~ 1.5$\sigma$ away from the retrieved dayside $T_B$.

Although these values of $A_g$ translate into the dayside MASCARA-1 reflecting a small fraction of the radiation incident upon it, this is a trend seen across the vast majority of hot Jupiters for which equivalent observations exist. Depending on which derivation of $A_g$ that we adopt (light curve fit, HELIOS retrieval or blackbody retrieval), the significance of the detection of reflected light ranges from $2-2.5\sigma$. We select the $A_g$ of $0.171_{-0.068}^{+0.066}$ derived from the light curve fit as our headline figure for this study, as the thermal and reflected components of flux are derived self-consistently using information from the full light curves, rather than the occultation depths alone. Wong et al. (2020a) and Wong et al. (2021) report that an observed positive correlation between $A_g$ and $T_d$ breaks down at $T_d \sim 3000$ K. Although the detection of reflected light in this study is relatively marginal, it is significantly higher than the estimated $A_g$ for other hot Jupiters with $T_d > 3000$ K such as WASP-18 b (Maxted et al. 2013, $A_g < 0.03$) and WASP-33 b (Zhang et al. 2018; von Essen et al. 2020, $A_g < 0.08$)<fnref n="18" /> This may indicate that the MASCARA-1 b dayside is not completely dominated by H⁻ absorption (Arcangeli et al. 2018; Lothringer et al. 2018), and that some highly reflecting species may not have dissociated close to the terminator.

## 8. Conclusions

We have presented a joint analysis of *Spitzer* and CHEOPS light curves of MASCARA-1 b, which has yielded the most precise constraints to date on a range of parameters including the period, epoch, eccentricity, stellar density, planetary radius, and orbital separation. The light curves of planets transiting fast rotating stars simultaneously encode information about the stellar inclination $i_\star$, the planetary inclination $i_p$, and the sky-projected spin-orbit angle $\lambda$, which in turn can be used to derive the true spin-orbit angle $\Psi$. However, the large degeneracies that exist between many of the transit parameters makes their accurate determination challenging. Doppler tomography observations do not resolve $\Psi$, but do generally allow $v \sin i_\star$ and $\lambda$ to be robustly measured. The use of priors informed by the tomography when fitting the transit light curves is a powerful tool to facilitate an accurate measurement of $\Psi$. Masuda (2015) recommend the placement of a Gaussian prior on both $v \sin i_\star$ and $\lambda$ from tomography when fitting such transits. As tomography also enables accurate measurements of $b$ that are less subject to degeneracies than the transits, we go one step further and suggest that Gaussian priors should be placed on $b$ in a similar way.

Our strategy of fitting both the CHEOPS and *Spitzer* transits with a common set of transit parameters also greatly reduced the effect of these degeneracies, measuring values of $v \sin i_\star$ and $\lambda$ consistent with the Doppler tomography even without the use of Gaussian priors, and allowed an independent measurement of $\Psi = 70.8_{-12.1}^{+10.8}$ deg. We propose that the acquisition of high precision transit light curves of fast rotators at multiple wavelengths may routinely provide a more robust determination of the transit parameters where Doppler tomography is not available, though this should be tested with upcoming observations. The use of CHEOPS to target hot Jupiters transiting fast rotating stars that have already been observed by TESS will provide a way to test whether observation in two overlapping optical bandpasses produces the same effect. As previous studies do not search for asymmetry in the *Spitzer* transits of fast rotators such as KELT-9 (Mansfield et al. 2020), Kepler-13A b (Shporer et al. 2014), and WASP-33 (Zhang et al. 2018), we suggest that the combined reanalysis of these light curves in conjunction with their corresponding TESS light curves may facilitate similar studies to the one we have performed for MASCARA-1.

When applying Gaussian priors on $b$, $v \sin i_\star$, and $\lambda$ informed by the Doppler tomography, we measured a more precise value of $\Psi = 72.1_{-2.4}^{+2.5}$ deg, which is the value we report. This result follows the trend of observed polar or near-polar orbits planets orbiting stars with $T_{\text{eff}} > 7000$ K, though measurements for a larger number of systems are needed to properly test theories describing the underlying physics.

In Sect. 6, we reported that the transit light curves show no evidence for significant orbital precession. As a result, we were not able to constrain the stellar gravitational quadrupole moment $J_2$: the term primarily dictating the rate of precession. The combination of the observational strategy and photometric precision afforded by the *Kepler* mission made it uniquely well-suited to search for orbital precession in hot Jupiters orbiting oblate stars—unfortunately the only target in the *Kepler* field for which this was possible was Kepler-13A b. The ability of CHEOPS to replicate such a study is also hampered by the gaps in the light curve, often limiting the ability of the impact parameter $b$ to be accurately determined from a single transit.

A strategy that future CHEOPS observations could adopt to mitigate against these issues could be to observe multiple transits of a suitable target across numerous years, ensuring that observations in a given observing season collectively sample all phases of the transit. When phase-folding all of the transits for each observing season, this strategy would facilitate search for orbital precession from year to year. Adopting this strategy for a sector of TESS observations allows Szabó et al. (2020) to confirm the continuing decrease in $b$ well after the *Kepler* mission had concluded. Although it would not come close to matching the time resolution achieved by *Kepler*, CHEOPS could observe optimal targets across a far wider area of the sky. Alongside MASCARA-1 b, this would include optimal targets such as KELT-9 b and WASP-189 b where the signatures of precession are expected to be larger and more easily detectable than for Kepler-13A b. TESS observations will supplement these searches, and for MASCARA-1 they are provisionally scheduled to occur in sector 52 (August 2022).

We also studied the atmosphere of MASCARA-1 b by measuring the occultation depths and phase-dependent flux modulation for both CHEOPS and *Spitzer*. Our light curve fitting suggests that MASCARA-1 b is within a small sample of hot Jupiters with an integrated dayside temperature $T_d > 3000$ K. The corresponding nightside temperature of $1720 \pm 330$ K and the derived heat redistribution factor of $0.571_{-0.078}^{+0.054}$ suggests

<footnote n="18">¹⁸ The upper limits are reported at $2\sigma$ by Wong et al. (2021).</footnote>
inefficient but existent recirculation of heat from the dayside. Our retrieval confirmed the assumption of our light curve model that the dayside emission spectrum is well described by a blackbody. While low, the measurement of $A_g = 0.171_{-0.068}^{+0.066}$ from our light curve model stands out compared to other ultra-hot Jupiters where $H^-$ absorption is expected to dominate in the dayside. The fact that our derived value of $A_s$ was somewhat larger than that of $A_B$ may suggest that MASCARA-1 b is less reflective at ultraviolet and near-infrared wavelengths than in the optical. The precision achieved for the CHEOPS occultations suggest that future CHEOPS observations of the MASCARA-1 b full-phase curve should significantly detect the planetary signal, which would allow the reflective properties at both optical and IR wavelengths to be separately determined.

To model the phase-dependent flux variations in the light curve, we used a novel, physically motivated phase curve parameterisation presented in M21 that generates a 2D temperature map to model the thermal emission and can describe any reflection law (Heng et al. 2021). This novel approach enabled us to separately derive the $A_g$, spherical albedo $A_s$, Bond albedo $A_B$, and phase integral $q$ of MASCARA-1. This approach to fitting multi-wavelength occultation and phase curve data in the optical and IR also provides a template for the analysis of spectrally resolved James Webb Space Telescope (JWST) data. Due to the limited phase coverage of the CHEOPS data, we fitted the data from both bandpasses with a single set of phase curve parameters, which implicitly assumed that they share a common albedo, and that the spectral energy distribution of the planet is well represented by a blackbody. Although a similar strategy may still be necessary when analysing occultation data, the use of a unique set of parameters per wavelength bin in the exquisite-precision spectrally resolved phase curves that JWST will acquire will allow for variations in thermal emission and reflection to be taken into account in the light curve fitting.

**Acknowledgements.** CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission. M.J.H. and Y.A. acknowledge the support of the Swiss National Fund under grant 200020_172746. S.H. gratefully acknowledges CNES funding through the grant 837319. D.K. acknowledges partial financial support from the Center for Space and Habitability (CSH), the PlanetS National Center of Competence in Research (NCCR), and the Swiss National Science Foundation and the Swiss-based MERAC Foundation. A.C.C. and T.G.W. acknowledge support from STFC consolidated grant number ST/M001296/1. P.M. acknowledges support from STFC research grant number ST/M001040/1. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2-190080). S.S. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 833925, project STAREX). A.Br. was supported by the SNSA. This work was supported by 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: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 and POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 and POCI-01-0145-FEDER-028953, PTDC/FIS-AST/28987/2017 and POCI-01-0145-FEDER-028987. O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia Maria de Maeztu-Centro de Astrobiología a (INTA-CSIC), as well as the support of the Generalitat de Catalunya/CERCA programme. The MOC activities have been supported by the ESA contract No. 4000124370. S.C.C.B. acknowledges support from FCT

through FCT contracts nr. IF/01312/2014/CP1215/CT0004. X.B., S.C., D.G., M.F., and J.L. acknowledge their roles as ESA-appointed CHEOPS science team members. This project was supported by the CNES. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (project FOUR ACES, grant agreement no. 724427). D.E. acknowledges financial support from the Swiss National Science Foundation for project 200021_200726. C.M.P. and M.F. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19, 174/18). D.G. gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 'Gaseousor rocky? Unveiling the nature of small worlds'. M.G. is an F.R.S.-FNRS Senior Research Associate. KGI is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme. She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection for the programme. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. Acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B-C33, as well as the support of the Generalitat de Catalunya/CERCA programme. S.G.S. acknowledges support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). This project has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) grants GINOP-2.3.2.15-2016-00003, K-119517, K-125015, and the City of Szombathely under Agreement No. 67.177-21/2016. V.V.G. is an F.R.S-FNRS Research Associate. We thank Taylor J. Bell for supplying and discussing the analysis of the Spitzer phase curve presented in Bell et al. (2021).

## References

Ahlers, J. P., Seubert, S. A., & Barnes, J. W. 2014, ApJ, 786, 131

Ahlers, J. P., Barnes, J. W., & Barnes, R. 2015, ApJ, 814, 67

Ahlers, J. P., Barnes, J. W., & Myers, S. A. 2019, AJ, 158, 88

Ahlers, J. P., Johnson, M. C., Stassun, K. G., et al. 2020a, AJ, 160, 4

Ahlers, J. P., Kruse, E., Colón, K. D., et al. 2020b, ApJ, 888, 63

Albrecht, S. H., Marcussen, M. L., Winn, J. N., Dawson, R. I., & Knudstrup, E. 2021, ApJ, 916, L1

Anderson, D. R., Temple, L. Y., Nielsen, L. D., et al. 2018, MNRAS, submitted [arXiv:1809.04897]

Angerhausen, D., DeLarme, E., & Morse, J. A. 2015, PASP, 127, 1113

Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30

Barnes, J. W. 2009, ApJ, 705, 683

Barnes, J. W., Linscott, E., & Shporer, A. 2011, ApJS, 197, 10

Barnes, J. W., Ahlers, J. P., Seubert, S. A., & Relles, H. M. 2015, ApJ, 808, L38

Bell, T. J., Dang, L., Cowan, N. B., et al. 2021, MNRAS, 504, 3316

Benomar, O., Masuda, K., Shibahashi, H., & Suto, Y. 2014, PASJ, 66, 94

Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109

Blackwell, D. E., & Shallis, M. J. 1977, MNRAS, 180, 177

Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18

Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5

Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021, A&A, 646, A157

Borsa, F., Lanza, A. F., Raspantini, I., et al. 2021, A&A, 653, A104

Borsato, L., Piotto, G., Gandolfi, D., et al. 2021, MNRAS, 506, 3810

Bourrier, V., Lovis, C., Beust, H., et al. 2018, Nature, 553, 477

Brothwell, R. D., Watson, C. A., Hébrard, G., et al. 2014, MNRAS, 440, 3392

Casasayas-Barris, N., Palle, E., Stangret, M., et al. 2021, A&A, 647, A26

Castelli, F., & Kurucz, R. L. 2003, in IAU Symposium, 210, Modelling of Stellar Atmospheres, eds. N. Piskunov, W. W. Weiss, & D. F. Gray, A20

Cegla, H. M., Lovis, C., Bourrier, V., et al. 2016, A&A, 588, A127

Chaplin, W. J., Sanchis-Ojeda, R., Campante, T. L., et al. 2013, ApJ, 766, 101

Chen, G., van Boekel, R., Wang, H., et al. 2014, A&A, 567, A8

Claret, A. 2016, A&A, 588, A15

Claret, A. 2021, Res. Notes Am. Astron. Soc., 5, 13

Claret, A., & Bloemen, S. 2011, A&A, 529, A75

Collier Cameron, A., Guenther, E., Smalley, B., et al. 2010, MNRAS, 407, 507

Cowan, N. B., & Agol, E. 2011, ApJ, 726, 82

Cizsmadia, S. 2020, MNRAS, 496, 4442

Dalal, S., Hébrard, G., Lecavelier des Étangs, A., et al. 2019, A&A, 631, A28

Deline, A., Hooton, M. J., Lendl, M., et al. 2021, A&A, in press [arXiv:2201.04518]

Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., 5, 775

Demory, B.-O., de Wit, J., Lewis, N., et al. 2013, ApJ, 776, L25
Demory, B.-O., Gillon, M., de Wit, J., et al. 2016, Nature, **532**, 207

Dorval, P., Talens, G. J. J., Otten, G. P. P. L., et al. 2020, A&A, **635**, A60

Dupret, M. A. 2001, A&A, **366**, 166

Dyudina, U., Zhang, X., Li, L., et al. 2016, ApJ, **822**, 76

Esteves, L. J., De Mooij, E. J. W., & Jayawardhana, R. 2015, ApJ, **804**, 150

Evans, T. M., Pont, F., Sing, D. K., et al. 2013, ApJ, **772**, L16

Foreman-Mackey, D. 2016, J. Open Source Softw., **1**, 24

Foreman-Mackey, D. 2018, Res. Notes Am. Astron. Soc., **2**, 31

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, **125**, 306

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, **154**, 220

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, **649**, A1

Gaudi, B. S., Stassun, K. G., Collins, K. A., et al. 2017, Nature, **546**, 514

Gibson, N. P. 2014, MNRAS, **445**, 3401

Gibson, N. P., Aigrain, S., Roberts, S., et al. 2012, MNRAS, **419**, 2683

Heng, K., & Li, L. 2021, ApJ, **909**, L20

Heng, K., & Workman, J. 2014, ApJS, **213**, 27

Heng, K., Morris, B. M., & Kitzmann, D. 2021, Nat. Astron., **5**, 1001

Henyey, L. G., & Greenstein, J. L. 1941, ApJ, **93**, 70

Herman, M. K., de Mooij, E. J. W., Huang, C. X., & Jayawardhana, R. 2018, AJ, **155**, 13

Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2019, Stat. Comput., **29**, 891

Hirano, T., Gaidos, E., Winn, J. N., et al. 2020, ApJ, **890**, L27

Hooton, M. J., Watson, C. A., de Mooij, E. J. W., Gibson, N. P., & Kitzmann, D. 2018, ApJ, **869**, L25

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

Howarth, I. D., & Morello, G. 2017, MNRAS, **470**, 932

Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, **635**, A24

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, **553**, A6

Johnson, M. C., Cochran, W. D., Albrecht, S., et al. 2014, ApJ, **790**, 30

Kamiaka, S., Benomar, O., & Suto, Y. 2018, MNRAS, **479**, 391

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., **90**, 773

Keating, D., Cowan, N. B., & Dang, L. 2019, Nat. Astron., **3**, 1092

Kipping, D. M. 2013, MNRAS, **435**, 2152

Kurucz, R. L. 2013, ATLAS12: Opacity sampling model atmosphere program

Leleu, A., Alibert, Y., Hara, N. C., et al. 2021, A&A, **649**, A26

Lendl, M., Csizmadia, S., Deline, A., et al. 2020, A&A, **643**, A94

Li, L., Jiang, X., West, R. A., et al. 2018, Nat. Commun., **9**, 3709

Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, **649**, A4

Lothringer, J. D., Barman, T., & Koskinen, T. 2018, ApJ, **866**, 27

Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars

Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, AJ, **153**, 56

Malik, M., Kitzmann, D., Mendonça, J. M., et al. 2019, AJ, **157**, 170

Mandel, K., & Agol, E. 2002, ApJ, **580**, L171

Mansfield, M., Bean, J. L., Stevenson, K. B., et al. 2020, ApJ, **888**, L15

Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, **835**, 77

Marley, M. S., Gelino, C., Stephens, D., Lunine, J. I., & Freedman, R. 1999, ApJ, **513**, 879

Masuda, K. 2015, ApJ, **805**, 28

Maxted, P. F. L., Anderson, D. R., Doyle, A. P., et al. 2013, MNRAS, **428**, 2645

Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2021, MNRAS in press, [arXiv:2111.08828]

May, E. M., & Stevenson, K. B. 2020, AJ, **160**, 140

McDonald, I., Zijlstra, A. A., & Boyer, M. L. 2012, MNRAS, **427**, 343

McLaughlin, D. B. 1924, ApJ, **60**, 22

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

Morris, S. L., & Naftilan, S. A. 1993, ApJ, **419**, 344

Morris, B. M., Delrez, L., Brandeker, A., et al. 2021, A&A, **653**, A173

Morris, B. M., Heng, K., Jones, K., et al. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202142135

Owens, N., de Mooij, E. J. W., Watson, C. A., & Hooton, M. J. 2021, MNRAS, **503**, L38

Pamyatnykh, A. A. 1999, Acta Astron., **49**, 119

Parviainen, H. 2015, MNRAS, **450**, 3233

Pearl, J. C., & Conrath, B. J. 1991, J. Geophys. Res., **96**, 18921

Perna, R., Heng, K., & Pont, F. 2012, ApJ, **751**, 59

Piskunov, N., & Valenti, J. A. 2017, A&A, **597**, A16

Queloz, D., Eggenberger, A., Mayor, M., et al. 2000, A&A, **359**, L13

Ramírez-Agudelo, O. H., Simón-Díaz, S., Sana, H., et al. 2013, A&A, **560**, A29

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telescopes Instrum. Syst., **1**, 014003

Rossiter, R. A. 1924, ApJ, **60**, 15

Russell, H. N. 1916, ApJ, **43**, 173

Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr., **90**, 054005

Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, Ap&SS, **316**, 83

Showman, A. P., & Polvani, L. M. 2011, ApJ, **738**, 71

Shporer, A., O'Rourke, J. G., Knutson, H. A., et al. 2014, ApJ, **788**, 92

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

Skilling, J. 2004, in Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, eds. R. Fischer, R. Preuss, & U. V. Toussaint, American Institute of Physics Conference Series, **735**, 395

Skilling, J. 2006, Bayesian Anal., **1**, 833

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, **131**, 1163

Southworth, J. 2011, MNRAS, **417**, 2166

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

Stefansson, G., Cañas, C., Wisniewski, J., et al. 2020, AJ, **159**, 100

Swayne, M. I., Maxted, P. F. L., Triaud, A. H. M. J., et al. 2021, MNRAS, **506**, 306

Szabó, G. M., Szabó, R., Benkő, J. M., et al. 2011, ApJ, **736**, L4

Szabó, G. M., Pál, A., Derekas, A., et al. 2012, MNRAS, **421**, L122

Szabó, G. M., Pribulla, T., Pál, A., et al. 2020, MNRAS, **492**, L17

Szabó, G. M., Gandolfi, D., Brandeker, A., et al. 2021, A&A, **654**, A159

Talens, G. J. J., Albrecht, S., Spronck, J. F. P., et al. 2017, A&A, **606**, A73

Talens, G. J. J., Albrecht, S., Spronck, J. F. P., et al. 2018, A&A, **613**, C2

Tsal, S.-M., Dobbs-Dixon, I., & Gu, P.-G. 2014, ApJ, **793**, 141

Valenti, J. A., & Piskunov, N. 1996, A&AS, **118**, 595

Van Grootel, V., Pozuelos, F. J., Thuillier, A., et al. 2021, A&A, **650**, A205

von Essen, C., Mallonn, M., Borre, C. C., et al. 2020, A&A, **639**, A34

von Zeipel, H. 1924, MNRAS, **84**, 665

Watanabe, N., Narita, N., & Johnson, M. C. 2020, PASJ, **72**, 19

Winn, J. N., Holman, M. J., Bakos, G. Á., et al. 2007, AJ, **134**, 1707

Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010, ApJ, **718**, L145

Wong, I., Shporer, A., Daylan, T., et al. 2020a, AJ, **160**, 155

Wong, I., Shporer, A., Kitzmann, D., et al. 2020b, AJ, **160**, 88

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

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, **140**, 1868

Zhang, M., Knutson, H. A., Kataria, T., et al. 2018, AJ, **155**, 83

Zhou, G., & Huang, C. X. 2013, ApJ, **776**, L35

<footnote n="1">¹ Physikalisches Institut, University of Bern, Gesellsschaftstrasse 6, 3012 Bern, Switzerland
e-mail: matthew.hooton@unibe.ch</footnote>

<footnote n="2">² Aix Marseille Univ, CNRS, CNES, LAM, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France</footnote>

<footnote n="3">³ Center for Space and Habitability, Gesellsschaftstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="4">⁴ Institute of Planetary Research, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany</footnote>

<footnote n="5">⁵ Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK</footnote>

<footnote n="6">⁶ Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, Versoix, Switzerland</footnote>

<footnote n="7">⁷ Astrophysics Group, Keele University, Staffordshire, ST5 5BG, UK</footnote>

<footnote n="8">⁸ Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="9">⁹ Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK</footnote>

<footnote n="10">¹⁰ Instituto de Astrofísica de Canarias, 38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="11">¹¹ Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, A-8042 Graz, Austria</footnote>

<footnote n="12">¹² Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden</footnote>

<footnote n="13">¹³ Instituto de Astrofísica e Ciencias do Espaco, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="14">¹⁴ Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote n="15">¹⁵ Department of Space, Earth and Environment, Onsala Space Observatory, Chalmers University of Technology, 439 92 Onsala, Sweden</footnote>

<footnote n="16">¹⁶ INAF, Osservatorio Astrofisico di Catania, Via S. Sofia 78, 95123 Catania, Italy</footnote>

<footnote n="17">¹⁷ Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="18">¹⁸ Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s/n, 08193 Bellaterra, Spain</footnote>
19 Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

20 Admatis, 5. Kandó Kálmán Street, 3534 Miskolc, Hungary

21 Depto. de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC campus, 28692 Villanueva de la Cañada (Madrid), Spain

22 Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

23 Université de Paris, Institut de physique du globe de Paris, CNRS, 75005 Paris, France

24 Centre for Mathematical Sciences, Lund University, Box 118, 22100 Lund, Sweden

25 Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

26 Space sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

27 INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

28 Leiden Observatory, University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands

29 Department of Space, Earth and Environment, Onsala space observatory, Chalmers University of Technology, 439 92 Onsala, Sweden

30 Dipartimento di Fisica, Universita degli Studi di Torino, via Pietro Giuria 1, 10125, Torino, Italy

31 University of Vienna, Department of Astrophysics, Türkenschanzstrasse 17, 1180 Vienna, Austria

32 Science and Operations Department – Science Division (SCI-SC), Directorate of Science, European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201-AZ Noordwijk, The Netherlands

33 Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15–17, Hungary

34 ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

35 Sydney Institute for Astronomy, School of Physics A29, University of Sydney, NSW 2006, Australia

36 IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, 75014 Paris, France

37 Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis bvd. Arago, 75014 Paris, France

38 Department of Astrophysics, University of Vienna, Tuerkenschanzstrasse 17, 1180 Vienna, Austria

39 INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

40 Institute of Optical Sensor Systems, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany

41 Dipartimento di Fisica e Astronomia “Galileo Galilei”, Universita degli Studi di Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy

42 ESTEC, European Space Agency, 2201AZ Noordwijk, The Netherlands

43 Center for Astronomy and Astrophysics, Technical University Berlin, Hardenberstrasse 36, 10623 Berlin, Germany

44 Institut für Geologische Wissenschaften, Freie Universität Berlin, 12249 Berlin, Germany

45 ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, 9700 Szombathely, Szent Imre h. u. 112, Hungary

46 MTA-ELTE Exoplanet Research Group, 9700 Szombathely, Szent Imre h. u. 112, Hungary

47 Ingenieurbüro Ulmer, Im Technologiepark 1, 15236 Frankfurt/Oder, Germany

48 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
Appendix A: Table of noise parameters

<table><caption>**Table A.1.** Similar to Table 5, but displaying the baseline and noise parameters. All parameters are dimensionless.</caption><thead><tr><th>Parameter</th><th>Prior</th><th>Best fit</th></tr></thead><tbody><tr><td>White noise Spitzer</td><td>lnσ<sub>w,S</sub> U(-12, -6)</td><td>-10.94<sup>+0.66</sup><sub>-0.92</sub></td></tr><tr><td>White noise CHEOPS</td><td>lnσ<sub>w,C</sub> U(-12, -6)</td><td>-9.348<sup>+0.66</sup><sub>-0.072</sub></td></tr><tr><td>GP hyperparameter</td><td>lnS<sub>0</sub> U(-30, 0)</td><td>-23.0<sup>+1.6</sup><sub>-1.0</sub></td></tr><tr><td>GP hyperparameter</td><td>lnω<sub>0</sub> U(-15, 15)</td><td>1.88<sup>+0.31</sup><sub>-0.45</sub></td></tr><tr><td>GP hyperparameter</td><td>lnQ</td><td>fixed</td></tr><tr><td>Constant 1</td><td>c<sub>1</sub> U(0, 2)</td><td>1.00673 ± 0.00021</td></tr><tr><td>y centroid 1</td><td>y<sub>1</sub> U(-1, 1)</td><td>0.0155<sup>+0.0015</sup><sub>-0.0014</sub></td></tr><tr><td>x<sup>2</sup> centroid 1</td><td>x<sub>1</sub><sup>2</sup> U(-1, 1)</td><td>-0.00116<sup>+0.00049</sup><sub>-0.00048</sub></td></tr><tr><td>y<sup>2</sup> centroid 1</td><td>y<sub>1</sub><sup>2</sup> U(-1, 1)</td><td>0.0063<sup>+0.0031</sup><sub>-0.0028</sub></td></tr><tr><td>xy centroid 1</td><td>xy<sub>1</sub> U(-1, 1)</td><td>-0.0 ± 0.001</td></tr><tr><td>x FWHM 1</td><td>xw<sub>1</sub> U(-1, 1)</td><td>-0.00426<sup>+0.00069</sup><sub>-0.00065</sub></td></tr><tr><td>y centroid 1</td><td>y<sub>1</sub> U(-1, 1)</td><td>-0.0083<sup>+0.0012</sup><sub>-0.0013</sub></td></tr><tr><td>y<sup>2</sup> FWHM 1</td><td>y<sub>1</sub><sup>2</sup> U(-1, 1)</td><td>-0.00044 ± 0.00058</td></tr><tr><td>Constant 2</td><td>c<sub>2</sub> U(0, 2)</td><td>1.00105<sup>+0.00054</sup><sub>-0.00051</sub></td></tr><tr><td>x centroid 2</td><td>x<sub>2</sub> U(-1, 1)</td><td>-0.0027<sup>+0.00037</sup><sub>-0.00036</sub></td></tr><tr><td>y centroid 2</td><td>y<sub>2</sub> U(-1, 1)</td><td>0.01109<sup>+0.0009</sup><sub>-0.00086</sub></td></tr><tr><td>x FWHM 2</td><td>xw<sub>2</sub> U(-1, 1)</td><td>-0.00519<sup>+0.00063</sup><sub>-0.00056</sub></td></tr><tr><td>y centroid 2</td><td>y<sub>2</sub> U(-1, 1)</td><td>0.0101<sup>+0.0013</sup><sub>-0.0012</sub></td></tr><tr><td>xy FWHM 2</td><td>xwy<sub>2</sub> U(-1, 1)</td><td>-0.00221<sup>+0.00078</sup><sub>-0.00082</sub></td></tr><tr><td>Constant 3</td><td>c<sub>3</sub> U(0, 2)</td><td>0.99927<sup>+0.00039</sup><sub>-0.00041</sub></td></tr><tr><td>x centroid 3</td><td>x<sub>3</sub> U(-1, 1)</td><td>-0.00399<sup>+0.00028</sup><sub>-0.00029</sub></td></tr><tr><td>y centroid 3</td><td>y<sub>3</sub> U(-1, 1)</td><td>-0.00285<sup>+0.00043</sup><sub>-0.00044</sub></td></tr><tr><td>y<sup>2</sup> centroid 3</td><td>y<sub>3</sub><sup>2</sup> U(-1, 1)</td><td>0.0096<sup>+0.0013</sup><sub>-0.0012</sub></td></tr><tr><td>xy centroid 3</td><td>xy<sub>3</sub> U(-1, 1)</td><td>-0.00175<sup>+0.00083</sup><sub>-0.00089</sub></td></tr><tr><td>x FWHM 3</td><td>xw<sub>3</sub> U(-1, 1)</td><td>-0.00279<sup>+0.00081</sup><sub>-0.00079</sub></td></tr><tr><td>x<sup>2</sup> FWHM 3</td><td>xw<sub>3</sub><sup>2</sup> U(-1, 1)</td><td>0.00107<sup>+0.00083</sup><sub>-0.00087</sub></td></tr><tr><td>y<sup>2</sup> FWHM 3</td><td>y<sub>3</sub><sup>2</sup> U(-1, 1)</td><td>0.0002<sup>+0.00057</sup><sub>-0.00062</sub></td></tr><tr><td>xy FWHM 3</td><td>xwy<sub>3</sub> U(-1, 1)</td><td>-0.0007<sup>+0.0013</sup><sub>-0.0011</sub></td></tr><tr><td>Constant 4</td><td>c<sub>4</sub> U(0, 2)</td><td>1.000073<sup>+0.000028</sup><sub>-0.000055</sub></td></tr><tr><td>Glint 4</td><td>g<sub>4</sub> U(0, 2)</td><td>1.01 ± 0.12</td></tr><tr><td>Constant 5</td><td>c<sub>5</sub> U(0, 2)</td><td>0.999919<sup>+0.000048</sup><sub>-0.000053</sub></td></tr><tr><td>Glint 5</td><td>g<sub>5</sub> U(0, 2)</td><td>0.97<sup>+0.12</sup><sub>-0.13</sub></td></tr><tr><td>Constant 6</td><td>c<sub>6</sub> U(0, 2)</td><td>0.999938<sup>+0.000056</sup><sub>-0.00004</sub></td></tr><tr><td>Sky background 6</td><td>bg<sub>6</sub> U(-1, 1)</td><td>0.000415<sup>+0.000087</sup><sub>-0.000079</sub></td></tr><tr><td>Glint 6</td><td>g<sub>6</sub> U(0, 2)</td><td>1.05<sup>+0.18</sup><sub>-0.17</sub></td></tr><tr><td>Constant 7</td><td>c<sub>7</sub> U(0, 2)</td><td>1.000174<sup>+0.00011</sup><sub>-0.00005</sub></td></tr><tr><td>Glint 7</td><td>g<sub>7</sub> U(0, 2)</td><td>1.108<sup>+0.089</sup><sub>-0.088</sub></td></tr></tbody></table>
**Appendix B: Corner plots of asymmetric transit fits**

<figure><img src="image_14.png" /><figcaption>**Fig. B.1.** Corner plot of the *CHEOPS* transit-only fit without the placement of Gaussian priors on *b*, *v sin i*, and *λ* informed by the Doppler tomography. The strategy is described in full in Section 5 and the corresponding results are displayed in Table 3.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. B.2.** Corner plot of the CHEOPS+Spitzer transit-only fit without the placement of Gaussian priors on *b*, *v sin i*, and *λ* informed by the Doppler tomography. The strategy is described in full in Section 5 and the corresponding results are displayed in Table 3.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

**Fig. B.3.** Corner plot of the *CHEOPS* transit-only fit with the placement of Gaussian priors on *b*, *v sin i*, and *λ* informed by the Doppler tomography. The strategy is described in full in Section 5 and the corresponding results are displayed in Table 3.
<figure><img src="image_17.png" /><figcaption>**Fig. B.4.** Corner plot of the CHEOPS+Spitzer transit-only fit with the placement of Gaussian priors on *b*, *v sin i*, and *λ* informed by the Doppler tomography. The strategy is described in full in Section 5 and the corresponding results are displayed in Table 3. This approach of using both the CHEOPS and Spitzer light curves and using priors from tomography is the one adopted for the fit to the full set of light curves presented in Section 7.</figcaption></figure>