© ESO 2022

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⁻¹ (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 emitted from the surface of the star would vary with local effective 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 temperature 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 λ, the fact that the stellar inclination i_⋆ is degenerate with the stellar rotational velocity v means that a precise measurement of the true, three-dimensional spin–orbit angle Ψ 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_⋆. A common method of breaking the aforementioned degeneracy is the use of asteroseismology to estimate i_⋆ (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° < i_⋆ < 80°. 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_⋆, and λ. This facilitates a measurement of Ψ 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 λ 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¹ both comparable to that of Kepler, this strategy allows CHEOPS to extend photometric measurements of Ψ 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 λ 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 (v sin i_⋆ = 109 ± 4 km s⁻¹), and that its planetary companion orbits in a near-polar orbit (λ = 69.5 ± 3.0°). With Gaia G_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 analysedjointly 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_eff, average radius R_⋆, mass M_⋆, and age t_⋆.

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_eff from the broad line wings of Hα 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⁻¹ from T17 and obtained T_eff = 7490 ± 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 [Fe/H] = 0.15 ± 0.15.

We determined R_⋆ 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; Gaia G, G_BP, and G_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± 0.022R_⊙.

We inferred M_⋆ and t_⋆ 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_eff, [Fe/H], and R_⋆ 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 indetail in Bonfanti et al. (2021), we combined the probability distributions of both M_⋆ and t_⋆ 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 δ Scuti- and γ 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 δ Scuti 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². 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 combinationof 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 discriminatebetween 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 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 depth materialising due to an increase in the brightness from the MASCARA-1 b dayside.

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.

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

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 y2 detector position in the baseline model. The violet lines display models instead including time in the baseline model.

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 aspart of programme 14059 with PI Jacob Bean, which we downloaded from the Spitzer Heritage Archive³. 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 signals, we created a custom model to simultaneously fit the light curves based on the existing tools available in pycheops⁴ (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.

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.

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⁵ (Parviainen 2015), which implements the widely used model first presented in Barnes (2009)⁶. 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(\vartheta )}{T_{\text{pole}}}={\left(\frac{g(\vartheta )}{g_{\text{pole}}}\right)}^{\beta },$$(1)

where ϑ is the stellar latitude, g is local surface gravity, g is polar surface gravity, and β 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\text{}\left(\vartheta \right)=-\frac{G{M}_{\star }}{r_{\vartheta }^2}\widehat{u_r}+{\left(\frac{2\pi }{P_{\text{rot}}}\right)}^2{R}_{\star }\sin \text{}\left(\vartheta \right)\widehat{u_{\perp }}.$$(2)

Here, G is the gravitational constant, M_⋆ is the stellar mass, P_rot is the stellar rotational period, R_⋆ is the stellar equatorial radius, and $\widehat{u_r}$ and $\widehat{u_{\perp }}$ are unit vectors radially from the centre and perpendicular to the stellar spin axis, respectively. r_ϑ is the latitute-dependent stellar radius accounting for the stellar oblateness f_⋆, and is equal to R_⋆ at the equator and R_⋆(1 − f_⋆) at the poles. In turn, f_⋆ is given by $$f_{\star }=\frac{3\pi }{2G{\rho }_{\star }{P}_{\text{rot}}^2},$$(3)

where ρ_⋆ 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 computed by integrating the stellar flux as a function of time over the filter response function. With T(ϑ) 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⁷ (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₀,

• transit depth ${\delta }_{\text{trans}}={(R_{\text{p}}/R_{\star })}^2$,

• stellar density ρ_⋆,

• impact parameter b,

• eccentricity e and argument of peristron ω_peri, parameterised as ecosω_peri and esinω_peri,

• and the quadratic limb darkening parameters q₁ and q₂ (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 λ: 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 measureorbital plane anti-clockwise with respect to the stellar equatorial plane (sometimes termed α), as opposed to the reverse direction (sometimes termed β).

• 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 ϕ – the angle between the equatorial plane and the observer’s line of sight – by the equation i_⋆ = 90 deg − ϕ.

• stellar rotational velocity projected onto the observer’s line of sight v sin i_⋆,

• gravity darkening coefficient β,

• and stellar polar temperature T_pole.

This parameterisation differs from that of Masuda (2015) by the use of v sin i_⋆ instead of the stellar rotation frequency f_rot, and ρ_⋆ instead of M_⋆: all of which are related by the equation $$\frac{v\sin i_{\star }}{2\pi f_{\text{rot}}}={\left(\frac{3M_{\star }}{4\pi {\rho }_{\star }}\right)}^{\frac{1}{3}}\sin i_{\star }.$$(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 < λ < 90 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_⋆, λ} and {−b, −i_⋆, −λ}. By applying the constraint of b > 0 and letting i_⋆ vary in the full range of − 90 deg < i_⋆ < 90 deg⁸, 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 θ and longitude ϕ. This is given by the equation (1) of M21, $$T\left(\theta ,\varphi \right)=\overline{T}\left(1+{\sum }_{m,\mathcal{l}}^{{\mathcal{l}}_{\text{max}}}{h}_{m,\mathcal{l}}(\theta ,\varphi )\right),$$(5)

where $\overline{T}$ is the mean background temperature. The h_m,ℓ basis functions describe perturbations to the mean background temperature and are given by Eq. (258) of Heng & Workman (2014), $$\begin{array}{ll}{h}_{m,\mathcal{l}}=\hfill & \frac{C_{m,\mathcal{l}}}{{\omega }_{\text{drag}}^2{\alpha }^4+m^2}{e}^{-{\tilde{\mu }}^2/2}[\mu m{H}_{\mathcal{l}}\cos (m\varphi )\hfill \\ \hfill & +\alpha {\omega }_{\text{drag}}\left(2l{H}_{\mathcal{l}-1}-\tilde{\mu }{H}_{\mathcal{l}}\right)\sin (m\varphi )]\hfill \end{array}$$(6)

where α is a dimensionless fluid number that depends on the Rossby and Prandtl numbers and ω_drag is the dimensionless drag frequency. M21 demonstrated that synthetic temperature maps generated by general circulation models (GCMs) are well-fitted by α = 0.6 and ω_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,ℓ 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,ℓ. Other quantities include μ = cosθ and $\tilde{\mu }=\alpha \mu$. The physicist’s Hermite polynomials are represented by H_ℓ. M21 further showed that one needs to specify a hot spot offset Δϕ 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 ξ (normalised such that transits occur at ξ = ±π) according toCowan & Agol (2011): $$\frac{F_{\text{p}}}{F_{\star }}=\frac{1}{\pi I_{\star }}{\left(\frac{R_{\text{p}}}{R_{\star }}\right)}^2{\displaystyle {\int }_0^{\pi }{\displaystyle {\int }_{-\xi -\pi /2}^{-\xi +\pi /2}{I}_{\text{p}}}}(\theta ,\varphi )\cos (\varphi +\xi )\sin \text{{^}2(\theta )\text{}d}\varphi \text{d}\theta ,$$(7)

with planetary intensity I_p given by $${\displaystyle \int \lambda }{\mathcal{F}}_{\lambda }{\mathcal{B}}_{\lambda }(T(\theta ,\varphi ))\text{d}\lambda$$(8)

for a filter with photon flux response function ${\mathcal{F}}_{\lambda }$ in units of electrons per photon, where ${\mathcal{B}}_{\lambda }$ is a Planck function. For I_⋆, we used a PHOENIX spectrum interpolated according to the values of T_eff, log g, and [Fe/H] listed in Table 1.

Alongside other parameters already defined in Sect. 4.1, for ℓ_max = 1⁹ the thermal phase modulation can be modelled using three parameters: C_1,1 ¹⁰, Δϕ, and $\overline{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, Δϕ is a constant phase offset of the sinusoid, and $\overline{T}$ denotes the DC constant offset term from zero-flux. We allowed Δϕ and $\overline{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 ε, 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 Ψ_ph. These quantities are in turn functions of the single-scattering albedo ω_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 Ψ_ph for a homogeneously reflective sphere. Upon obtaining Ψ_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_{\text{1},\text{1}},\overline{T},\Delta \varphi ,{\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 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 lightcurves (Gibson 2014). The GP is parameterised by the frequency of the undamped oscillator ω₀, S₀ is proportional to the power in the power spectral density at ω = ω₀, 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 dynesty 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 $\mathcal{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 $\mathcal{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 β (not to be confused with the gravity darkening coefficient β) 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 β.

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 λ 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 λ = 58.6 ± 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¹²). 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 λ 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 Ψ, which as described in Sect. 1 is difficult to measure using other techniques. This is given by the relation: $$\cos \Psi =\cos i_{\star }\cos i_{\text{p}}+\sin i_{\star }\sin i_{\text{p}}\cos \lambda ,$$(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 λ. Thus, an accurate measurement of Ψ 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 ρ_⋆) and T_eff. We can infer the values of q₁ and q₂ from previous studies based on the stellar type and wavelength dependency (Claret & Bloemen 2011; Claret 2021). Similarly we can infer the value of β (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 vsini_⋆, and λ. 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 sini_⋆, and λ 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 sini_⋆ 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 λ (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 β. This is because as i_⋆ → 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 β to fit any asymmetry, and vice versa. Even when trying to mitigate against this by placing a Gaussian prior on β, the fits regularly returned results with β at an unphysically high value. For this reason, we fixed β = 0.199 according toClaret (2016) for later runs.

We also placed a Gaussian prior R_⋆ 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_pole = 7490 K according tothe 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 ρ_⋆ according to our determinations of R_⋆ and M_⋆, 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₁ and q₂ vary uniformly. Finally, we also let δ_trans and i_⋆ vary uniformly within wide limits, and P and t₀ 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 sini_⋆, and λ. 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. δ_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_⋆ 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σ 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 then. 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. 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 < λ < 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 < λ < 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_⋆ due to the symmetry of the star about the equatorial plane.

Secondly, in most circumstances Doppler tomography facilitates accurate determination of b, v sini_⋆, and λ, but not Ψ due to the lack of information about i_⋆. 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_⋆, which enables an accurate derivation of Ψ. 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 β 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{{\mathcal{Z}}_{\text{Barnes}}}{{\mathcal{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σ 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 β 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₂ 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 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₁₄) and λ. Previous measurements of J₂ 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⁻⁴. Adopting the same value for MASCARA-1, we would anticipate observing d b∕dt ~ 0.002 and d λ∕dt ~ 0.34 deg yr⁻¹. The precision of the data will prevent us to detect variations at this low level, but J₂ 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 λ to median values in the CHEOPS + Spitzer fit with prior informed by the Doppler tomography as described in the previous section, and let b, λ 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 λ 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 sini_⋆, λ to the values measured by the Doppler tomography in T17, and let P, t₀, δ_trans, and ρ_⋆ and the limb darkening parameters per instrument (q₁ and q₂) 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σ of the measurement of b = 0.122 ± 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σ 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 Spitzer and second CHEOPS transit. We discuss strategies for a possible future study to explore this in greater depth in Sect. 8.

Fig. 5 Phase-folded Spitzer and CHEOPS photometry centred on the transits and occultations, using the same format as the previouslight 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.

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_⋆ and λ 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⁻¹ and $\lambda =-{69.2}_{-3.4}^{+3.1}$ deg are both within 2σ of the Doppler tomography¹⁴. Unlike the Doppler tomography, the transit light curve encodes information about the stellar inclination i_⋆, which we measured to be ${55.5}_{-2.9}^{+2.3}$ deg. As this breaks the degeneracy between v and i_⋆, we measured the stellar rotation period P_rot to be ${0.853}_{-0.016}^{+0.017}$ days. This rapid rotation causes the stellar equator to bulge to a radius 104.39 ± 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 λ (top panel) and Ψ (bottom panel) are shown in Fig. 7 in the context of similar measurements collected in the TEPCat orbital obliquity catalogue¹⁵ (Southworth 2011). Spectroscopic methods have facilitated many more measurements of λ than Ψ, with stars where T_eff < 7000 K most likely to occupy aligned, prograde orbits. Stars with T_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_eff comparable to measurements for Kepler-13A b 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.

Recently, Albrecht et al. (2021) note that the majority of planets with misaligned orbits occupy nearly polar orbits in the range Ψ = 80–120 deg. Like Kepler-13A b, while relatively close to polar the Ψ of MASCARA-1 b falls outside this particularly clustered range. Although not as massive as Kepler-13A b, with M_p = 3.7 ± 0.9M_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_eff > 7000 K will be necessary to reveal more detailed trends.

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 λ. The {b, λ} 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, λ} 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, λ} 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⋆, λ} and another with {−b, −i⋆, −λ}. The {b, λ} pairs with positive i⋆ (north pole pointing towards the observer) are shown in the left panel, and the {b, λ} 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 < λ < 90 deg, and corresponding solutions are found for all of the transit fits presented in this paper.

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,ℓ basis functions, which replicate the shape of the global temperature distributions produced by many GCMs (e.g. Showman & Polvani 2011). As ω_scat is fixed to 4.5 in the fit, the hotspot offset Δϕ is roughly equivalent to the offset between the phase of superior conjunction and maximum flux, which with $\Delta \varphi ={2.0}_{-9.4}^{+8.9}$ deg is consistentwith zero.

A corner plot showing the correlations between the phase curve parameters is shown in Fig. 9. The negative correlation between lnc_1,1 and $\overline{T}$ is analogousto 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 ω_scat, which scalesthe 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_{\text{g}}={0.171}_{-0.068}^{+0.066}$ constitutes a 2.5σ 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.

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σ higher than the Bell et al. (2021) equivalents. Although different approaches were 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.

Fig. 7 All published measurements of the absolute value of the sky projected spin–orbit angle |λ| (top panel) and true spin–orbit angle Ψ (bottom panel) as a function of Teff, according to TEPCat. All planets with measurements of Ψ with Teff > 7000 K are marked specifically. In cases where more than one measurement exists for a given planet, they are all shown. Some values of λ were converted to fall in the range − 180 < λ < 180 deg.

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.

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).

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.

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 nightside 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_{\text{B}}=1-{\left(\frac{a}{R_{\star }}\right)}^2\frac{{\displaystyle {\int }_{-\pi }^{\pi }{\displaystyle {\int }_0^{\pi }{F}_{\text{p}}}}\sin \theta \text{d}\theta \text{d}\varphi }{\pi {\sigma }_{\text{SB}}{T}_{\text{eff}}^4},$$(10)

where σ_SB is the Stefan-Boltzmann constant. The redistribution factor ε is simply the ratio of the flux from the integrated nightside and dayside (Perna et al. 2012), $$\epsilon =\frac{{\displaystyle {\int }_{\pi /2}^{\pi }{\displaystyle {\int }_0^{\pi }{F}_{\text{p}}}}\sin \theta \text{d}\theta \text{d}\varphi +{\displaystyle {\int }_{-\pi }^{-\pi /2}{\displaystyle {\int }_0^{\pi }{F}_{\text{p}}}}\sin \theta \text{d}\theta \text{d}\varphi }{{\displaystyle {\int }_{-\pi /2}^{\pi /2}{\displaystyle {\int }_0^{\pi }{F}_{\text{p}}}}\sin \theta \text{d}\theta \text{d}\varphi }.$$(11)

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

The derived values of ε and A_B using both the 0D approach from Cowan & Agol (2011) and 2D approach from M21 are displayed 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_B, which suggests that it could be underestimated due to spectral features in the emission spectrum of MASCARA-1 b within the Spitzer/IRAC 4.5 μm bandpass causing 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_⋆, integrated over all wavelengths (Marley et al. 1999), $$A_{\text{B}}=\frac{{\displaystyle {\int }_0^{\infty }{A}_{\text{s}}}\left(\lambda \right) F_{\star } \text{d}\lambda }{{\displaystyle {\int }_0^{\infty }{F}_{\star }} \text{d}\lambda }.$$(12)

In this study, each quantity is derived separately: A_B as describedin the previous paragraphs using the thermal component of the light curve model and A_s as the product of A_g and the phase integral 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 intrinsic properties of the atmosphere that do not formally involve the star. In practice, spherical and geometric albedos are commonly 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¹⁶ measurement of $A_{\text{s}}={0.266}_{-0.100}^{+0.097}$ is larger than our measurement of A_B from the 2Dtemperature 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 optical 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 bandpasses (i.e. A_s,out = 0), and assuming that optical wavelengths and longer are well-described by the A_s derived in our light curve model (i.e. $A_{\text{s},\text{in}}={0.266}_{-0.100}^{+0.097}$), using Eq. (12) we computed $A_{\text{B}}={0.184}_{-0.069}^{+0.067}$. Although still over double the A_B 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 temperatureof $T_{\text{d}}={3062}_{-68}^{+66}$ K confirms the position of the MASCARA-1 dayside amongst the hottest of any known exoplanet. We report a 3σ detection offlux from the nightside and a corresponding nightside integrated temperature of T_n = 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 nightside 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_n 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 analysis of the thermal and reflective properties, for the sake of computational efficiency it necessarily included some major simplifications. 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 atmosphere, we generated a grid of emission spectra using HELIOS¹⁷ (Malik et al. 2017, 2019) to jointly interpret the CHEOPS and Spitzer occultation depths.

We used a planetary radius of 1.597 R_J and an orbital distance of 0.040352 au derived from our light curve fitting, and 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 ε (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 onlythe 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).

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 ε, we parameterised them with a spline function for each filter.

We then used an MC algorithm to find the distribution of ε 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_⋆ 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 ε = 0.307 ± 0.088 (f = 0.538 ± 0.038). The spectra for the median value of ε 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 ± 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 ± 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 ± 53 K, which resulted in a predicted CHEOPS occultation depth of 90 ppm and an A_g of 0.127 ± 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_{\text{g}}={0.171}_{-0.068}^{+0.066}$ is in good agreement with both the retrieved HELIOS and blackbody values, and $T_{\text{d}}={3062}_{-68}^{+66}$ K ~ 1.5σ 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σ. 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 ~ 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)¹⁸ 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.

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.

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

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_⋆, the planetary inclination i_p, and the sky-projected spin–orbit angle λ, which in turn can be used to derive the true spin–orbit angle Ψ. However, the large degeneracies that exist between many of the transit parameters makes their accurate determination challenging. Doppler tomography observations do not resolve Ψ, but do generally allow vsini_⋆ and λ 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 Ψ. Masuda (2015) recommend the placement of a Gaussian prior on both vsini_⋆ and λ 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 sini_⋆ and λ 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_⋆, and λ 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_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₂: 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 toyear. 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 tooccur 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 ± 330 K and the derived heat redistribution factor of ${0.571}_{-0.078}^{+0.054}$ suggests 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_{\text{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 Internacionalizacã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 (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 projecthas 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 GrantNo. 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).

Appendix A Table of noise parameters

Appendix B Corner plots of asymmetric transit fits

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.

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.

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 tomo- graphy. The strategy is described in full in Section 5 and the corresponding results are displayed in Table 3.

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.

References

All Tables

All Figures

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.

In the text

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

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 y2 detector position in the baseline model. The violet lines display models instead including time in the baseline model.

In the text

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.

In the text

	Fig. 5 Phase-folded Spitzer and CHEOPS photometry centred on the transits and occultations, using the same format as the previouslight 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.
In the text

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 λ. The {b, λ} 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, λ} 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, λ} 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⋆, λ} and another with {−b, −i⋆, −λ}. The {b, λ} pairs with positive i⋆ (north pole pointing towards the observer) are shown in the left panel, and the {b, λ} 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 < λ < 90 deg, and corresponding solutions are found for all of the transit fits presented in this paper.

In the text

Fig. 7 All published measurements of the absolute value of the sky projected spin–orbit angle |λ| (top panel) and true spin–orbit angle Ψ (bottom panel) as a function of Teff, according to TEPCat. All planets with measurements of Ψ with Teff > 7000 K are marked specifically. In cases where more than one measurement exists for a given planet, they are all shown. Some values of λ were converted to fall in the range − 180 < λ < 180 deg.

In the text

	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.
In the text

	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).
In the text

	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.
In the text

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.

In the text

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

	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.
In the text

	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.
In the text

	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 tomo- graphy. The strategy is described in full in Section 5 and the corresponding results are displayed in Table 3.
In the text

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.

In the text