ABSTRACT

WASP-33b, a hot Jupiter around a hot star, is a rare system in which nodal precession has been discovered. We updated the model for the nodal precession of WASP-33b by adding new observational points. Consequently, we found a motion of the nodal precession spanning 11 yr. We present homogenous Doppler tomographic analyses of eight data sets, including two new data sets from TS23 and HIDES, obtained between 2008 and 2019, to illustrate the variations in the projected spin–orbit obliquity of WASP-33b and its impact parameter. We also present its impact parameters based on photometric transit observations captured by MuSCAT in 2017 and MuSCAT2 in 2018. We derived its real spin–orbit obliquity ψ, stellar spin inclination is, and stellar gravitational quadrupole moment J2 from the time variation models of the two orbital parameters. We obtained |$\psi = 108.19^{+0.95}_{-0.97}$| deg, |$i_\mathit{ s} = 58.3^{+4.6}_{-4.2}$| deg, and |$J_2=(1.36^{+0.15}_{-0.12}) \times 10^{-4}$|⁠. Our J2 value was slightly smaller than the theoretically predicted value, which may indicate that its actual stellar internal structure is different from the theoretical one. We derived the nodal precession speed |$\dot{\theta }=0.507^{+0.025}_{-0.022}$| deg yr−1, and its period |$P_{\mathrm{pre}}=709^{+33}_{-34}$| yr, and found that WASP-33b transits in front of WASP-33 for only ∼ 20 per cent of the entire nodal precession period.

1 INTRODUCTION

To date, WASP (Wide Angle Search for Planets; Collier Cameron et al. 2007), KELT (Kilodegree Extremely Little Telescope; Pepper et al. 2007), and TESS (Ricker et al. 2015) surveys have confirmed that there are 17 hot Jupiters around hot stars (>7000 K). Despite their small number, their orbital obliquities tend to have a wide range. This tendency indicates that they did not experience realignment by the tidal torque (Albrecht et al. 2012). Generally, hot host stars rotate rapidly, yielding oblateness larger than those of slowly rotating stars. A misaligned orbit and a fast-rotating star force faster precession. At present, there are only two planets whose nodal precessions have been detected: Kepler-13Ab (Szabó et al. 2012) and WASP-33b (Johnson et al. 2015). These planets are hot Jupiters, revolving in misaligned orbits around rapidly rotating hot stars.

WASP-33b, a hot Jupiter (Rp = 1.5RJ) around an A-type (Teff = 7430 ± 100 K) and rapidly rotating (Vsin is = 85.6 km s−1) star, is the only planet whose precession has been detected by Doppler tomography based on transit spectral data. This planet was discovered by Collier Cameron et al. (2010); its nodal precession has been detected by more than one Doppler tomographic measurement (Johnson et al. 2015; Iorio 2016; Watanabe, Narita & Johnson 2020; Borsa et al. 2021). The latest study on this topic (Borsa et al. 2021) derived the angle between the stellar spin axis and line of sight is = 90.11 ± 0.12 deg as well as the stellar gravitational quadrupole moment J2 = (6.73 ± 0.22) × 10−5 from the measurements; however, the value of is does not match that of Iorio (2016) within 3σ, whereas the value of J2 is also not consistent with that of Watanabe et al. (2020) within 3σ. Moreover, (Borsa et al. 2021) revealed that the change in the transit chord of WASP-33b is slightly overly complicated for hereafter estimating the position of the transit chord. Thus, we performed more additional observations and data sets to obtain more accurate values of is and J2, as well as a clearer forecasting of the transit chord.

In this paper, we report the nodal precession of WASP-33b for a longer period of 11 yr by considering additional Doppler tomographic measurements obtained through spectral transit and transit photometric observations. Transit photometry is an important method for measuring the impact parameter b and its change, as demonstrated by Szabó et al. (2012), who applied this to detect the nodal precession of Kepler-13Ab using Kepler’s photometric data sets. In Section 2, we present our spectral data sets for Doppler tomography and the steps for measuring the orbital parameters, i.e. b and the projected spin–orbit obliquity λ, the angle between the stellar spin axis and the planetary orbital axis. We explain how to handle our transit photometric data sets for determining b in Section 3. Moreover, we fit the values of these two parameters from the spectral and photometric measurements with the nodal precession model, which is described in Section 4. We display the behaviour of WASP-33b’s nodal precession and the derived parameters of WASP-33b and its host star in Section 5. In Section 6, we discuss the nodal precession results. Finally, we present our conclusions in Section 7.

2 DOPPLER TOMOGRAPHIC MEASUREMENT

We can simultaneously measure λ and b from the transit spectral data via Doppler tomography. When a planet passes in front of the stellar disc, a bump, referred to as a planetary shadow, appears in the stellar line profile. The orbital configuration of the planet can be derived from this shadow motion.

2.1 Observation data sets

We used eight spectroscopic data sets for WASP-33 around planetary transits. Table 1 lists the details of these transit spectral data sets.

One of them was obtained with the High Dispersion Spectrograph (HDS; Noguchi et al. 2002) at the 8.2 m Subaru telescope on 2011 October 19 ut. The other two data sets were obtained by the Harlan J. Smith Telescope (HJST) with the Robert G. Tull Coudé Spectrograph (TS23; Tull et al. 1995) at McDonald Observatory on 2008 November 12 ut and 2014 October 4 ut. These data sets were also used in our previous study Watanabe et al. (2020). To extract each line profile from each spectrum, we adopted least-squares deconvolution (LSD; Donati et al. 1997); the observed spectrum was regarded as a convolution of a line profile and a series of delta functions. Under this method, we obtained a list of the absorption lines from the Vienna Atomic Line Database (VALD; Kupka et al. 2000) to create a series of delta functions. We then derived each line profile and the error bars with the matrix calculations in Kochukhov, Makaganiuk & Piskunov (2010).

We then have analysed three spectroscopic data sets obtained on 2016 September 28 ut, 2018 January 12 ut, and 2019 January 2 ut from the high-resolution HARPS-N spectrograph (Cosentino et al. 2012), which is mounted at the Telescopio Nazionale Galileo (TNG). These three data sets have been extracted and published in Borsa et al. (2021), which were used in the extracted line profile series in this study.

Then, we also included two additional spectral data sets. One was an extracted data set taken by HJST/TS23 on 2016 December 11 ut. Data from 2016 also included 10 in-transit spectra. The other data set was obtained from the 188 cm telescope with HIgh Dispersion Echelle Spectrograph (HIDES; Izumiura 1999) at Okayama Astro-Complex (OAC) in Japan on 2019 December 27 ut. We utilized a wavelength range from 4980 Å to 6220 Å except for the Na D lines and wavelength regions around bad pixels. We reduced these spectral data by the subtracting bias and dark features, flat dividing, and performing 1D spectrum and wavelength calibration. We then utilized their continua, as well as the HDS process. We selected alp Leo, a rapidly rotating star (Vsin is ∼ 300 km s−1; Abt, Levato & Grosso 2002), to erase the Earth’s atmospheric absorption lines. We then shifted the spectra to the barycentric frame using pyraf (Science Software Branch at STScI 2012). Finally, we extracted each line profile of each exposure using the LSD.

2.2 Extracting planetary shadow

We derived a median line profile from all exposure data including in-transit for each epoch. We applied a median line profile, not a mean profile, because it avoided the effects of outliers. We subtracted the median line profile from each exposure line profile to calculate the time-series of the line profile residuals. Both the planetary shadow due to the transit of WASP-33b and a striped pattern due to the non-radial pulsations on the surface of WASP-33 (Collier Cameron et al. 2010) were present in the line profile residuals. Herrero et al. (2011) showed that the pulsation period was approximately 68 min based on photometric observations. However, determining the period from the Doppler tomographic results was difficult owing to the irregular patterns.

To extract only the planetary shadow, we applied a Fourier filtering technique (Johnson et al. 2015) because the planet is retrograde, while the pulsations are prograde. First, we performed a 2D Fourier transform. There were components derived from the pulsations in the Fourier space in the first and third quadrants, and the planetary shadow’s components in the second and fourth quadrants. Secondly, we created a filter in which we set unity in the two diagonal quadrants, including power from the planetary shadow, and zero in the other quadrants, including power from the pulsation with a Hann function between these quadrants. Finally, we multiplied the Fourier space by the filter and performed an inverse Fourier transform on the filtered Fourier space to extract the planetary shadow. Fig. 1 illustrates these procedures.

Doppler tomographic data sets and Fourier filters. Left column: observed residuals for the line profile series. The vertical dotted lines show v = 0 and ±vsin is. The bottom, middle, and upper horizontal dotted lines show the beginning, middle, and end of the WASP-33b’s transit, respectively. Middle column: Fourier spaces after Fourier transform of the residuals for the line profile series. These colour scales are shown as square roots. The faint narrow structure from the bottom right to the upper left is a component of WASP-33b’s planetary transit. In contrast, the bright-wide structure from the bottom left to the upper right is a pulsation component. Right column: filtered Fourier space such that only the transit component remains.
Figure 1.

Doppler tomographic data sets and Fourier filters. Left column: observed residuals for the line profile series. The vertical dotted lines show v = 0 and ±vsin is. The bottom, middle, and upper horizontal dotted lines show the beginning, middle, and end of the WASP-33b’s transit, respectively. Middle column: Fourier spaces after Fourier transform of the residuals for the line profile series. These colour scales are shown as square roots. The faint narrow structure from the bottom right to the upper left is a component of WASP-33b’s planetary transit. In contrast, the bright-wide structure from the bottom left to the upper right is a pulsation component. Right column: filtered Fourier space such that only the transit component remains.

2.3 Deriving parameters

To obtain the best-fitting values and uncertainties of the transit parameters, we adopted the Markov chain Monte Carlo (MCMC) method using the code emcee (Foreman-Mackey et al. 2013).

We modelled a planetary shadow via a convolution between a rotational broadening profile and a Gaussian line profile owing to intrinsic broadening, thermal broadening, and micro-turbulence. The detailed equations to derive the model of the planetary shadow are described in the appendix of Watanabe et al. (2020). We then applied the same filter to the model of the planetary shadow following the procedures described in Section 2.2.

We fitted the observed residuals of the five data sets to the models with 21 parameters using MCMC: the λ, b, and transit mid-time Tc of each epoch, Vsin is, Rp/Rs, a/Rs, two quadratic limb darkening coefficients, and the full width at half-maximum (FWHM) of the Gaussian line profile. The limb darkening coefficients were derived via the triangular sampling method reported in Kipping (2013), with q1, and q2. We estimated that q1 and q2 of HDS, TS23, and HIDES were equivalent. They were calculated from the stellar parameters, i.e. the effective temperature Teff, surface gravity log g, and metallicity. We set the priors of λ and b for all epochs and the FWHM as uniform functions while that of the priors of the other parameters were set as Gaussian priors. For the values and widths of the Gaussian priors, we set the priors for Rp/Rs and a/Rs based on the values and uncertainties reported in Kovács et al. (2013); the priors of each Tc of each epoch from Porb in von Essen et al. (2014) and T0 in von Essen et al. (2019); those of q1, and q2 calculated via PyLDTk (Husser et al. 2013; Parviainen & Aigrain 2015); and those of Vsin is from Johnson et al. (2015).

For the fitting, we maximized the logarithm of the posterior probability, ln Lpost:
(1)
where Oi is the data, Ci is the model, σi is the error for the ith data point, p is the parameter value at the gained iteration of the Markov chain, μ is the value from the literature, and s is the uncertainty from the literature. Index j denotes the parameters of the Gaussian priors. We set the range of the uniform prior of each λ to −180 deg < λ < −90 deg and that of each b as −1 < b < 1. To converge these parameter values, we ran 4000 steps, cut-off the first 2000 steps as burn-in, and iterated this set 100 times. Figs A1 and A2 in Appendix  A plot the posterior distributions.
Fitting of the filtered residual data by MCMC. Same type of colour scale as that in the left column of Fig. 1. Left column: residual data remaining as a planetary shadow. Middle column: filtered models of a planetary shadow using the best-fitting values. Right column: the difference between the first row and the second rows.
Figure 2.

Fitting of the filtered residual data by MCMC. Same type of colour scale as that in the left column of Fig. 1. Left column: residual data remaining as a planetary shadow. Middle column: filtered models of a planetary shadow using the best-fitting values. Right column: the difference between the first row and the second rows.

Light curves of MuSCAT in 2017 with the original cadences (grey points) and 1-min bins (coloured points). Top row: light curves from the photometric observation data. The black solid lines represent the models from the MCMC fitting. Middle row: light curves subtracted using the Gaussian process. Bottom row: residuals between the observed data and model data.
Figure 3.

Light curves of MuSCAT in 2017 with the original cadences (grey points) and 1-min bins (coloured points). Top row: light curves from the photometric observation data. The black solid lines represent the models from the MCMC fitting. Middle row: light curves subtracted using the Gaussian process. Bottom row: residuals between the observed data and model data.

Similar light curves as Fig. 3, but for the MuSCAT2 data set.
Figure 4.

Similar light curves as Fig. 3, but for the MuSCAT2 data set.

MCMC corner plots for the architecture angles, ψ, θ2008, and is, and stellar quadrupole moment J2 of the WASP-33b System.
Figure 5.

MCMC corner plots for the architecture angles, ψ, θ2008, and is, and stellar quadrupole moment J2 of the WASP-33b System.

Changes in λ (upper row) and b (lower row). The left and right columns show the short- and long-term changes, respectively. The x-axis of the right column is the time in years from the epoch of 2008. The blue circles show values from the HJST/TS23 data sets, the green triangles are values from Subaru/HDS, the red squares are values from OAC/HIDES, the cyan line is the value range from OAC/MuSCAT, the magenta line is the value range from TCS/MuSCAT2, and the black solid lines represent the model. In the bottom-right figure, the two black-dashed lines show the edges of the stellar disc of WASP-33b.
Figure 6.

Changes in λ (upper row) and b (lower row). The left and right columns show the short- and long-term changes, respectively. The x-axis of the right column is the time in years from the epoch of 2008. The blue circles show values from the HJST/TS23 data sets, the green triangles are values from Subaru/HDS, the red squares are values from OAC/HIDES, the cyan line is the value range from OAC/MuSCAT, the magenta line is the value range from TCS/MuSCAT2, and the black solid lines represent the model. In the bottom-right figure, the two black-dashed lines show the edges of the stellar disc of WASP-33b.

3 PHOTOMETRIC MEASUREMENT

3.1 Photometric observations of WASP-33b’s transit

We observed the transit of WASP-33b by photometry using two instruments with multicolour simultaneous cameras: Multicolor Simultaneous Camera for studying Atmospheres of Transiting exoplanets (MuSCAT; Narita et al. 2015) on the 188 cm telescope at OAC and MuSCAT2 (Narita et al. 2019) on the Telescopio Carlos Sánchez (TCS) 1.52 m telescope at the Teide Observatory (OT). MuSCAT has three channels for the |$g^{^{\prime }}_{2}$| (400–550 nm), |$r^{^{\prime }}_{2}$| (550–700 nm), and zs,2 (820–920 nm) bands. In contrast, MuSCAT2 contains four channels for the |$g^{^{\prime }}_{2}$|⁠, |$r^{^{\prime }}_{2}$|⁠, zs,2, and |$i^{^{\prime }}_{2}$| (700–820 nm) bands. These bands are the Astrodon Photometrics Generation 2-type Sloan filters. We obtained the data set from MuSCAT on 2017 November 5 ut. We set the exposure times of the |$g^{^{\prime }}_{2}$|⁠, |$r^{^{\prime }}_{2}$|⁠, and zs,2 bands for 4, 4, and 10 s, respectively. We also obtained the data set from MuSCAT2 on 2018 October 11 ut. The exposure time of the |$g^{^{\prime }}_{2}$|⁠, |$r^{^{\prime }}_{2}$|⁠, |$i^{^{\prime }}_{2}$|⁠, and zs,2 was 3, 2, 5, and 12 s, respectively.

To produce the light curves of WASP-33b, aperture photometry was performed using the pipeline proposed by Fukui et al. (2011). This process determines the stellar barycentre in every frame. We used BD+36 488, the second brightest star image in the frame, as a companion star to correct the atmospheric extinction. Next, the pipeline calculated the shift in the stellar position relative to the reference frame. This approach photometers the target star and a comparison star with a fixed aperture radius. Here, we set the aperture radii to 36 and 40 pixels for MuSCAT and MuSCAT2, respectively. After the sky background in the torus area centred on the stellar barycentre, the pipeline subtracted the sky background from WASP-33b’s flux and that of the comparison. Finally, the WASP-33b light curve was obtained by dividing its flux by the comparison star flux.

3.2 Light-curve fitting

To measure WASP-33b’s impact parameter b in 2017 and 2018, we constructed light-curve models with a Gaussian process using the python code exoplanet (Foreman-Mackey et al. 2020). WASP-33b’s light curve not only showed dimming by the transit, but also a short sinusoidal-wave-like feature due to the stellar pulsations. Thus, following a previous study by Johnson et al. (2015), we applied a Matern 3/2 kernel Kker, whose element is expressed as follows:
(2)
for the Gaussian process. i and j denote the orders of the photometric observation’s data, ti and tj are the observation times, α and l are the hyper parameters indicating the amplitude and time-scale of the stellar variations, respectively, and σ is the error of data point i.
We then fitted the light curves on two epochs to the models with the following 30 parameters using MCMC: baseline B for each light curve, b and Tc of each epoch, two quadratic limb darkening coefficients u1, u2, and α of each band, Rp/Rs, Porb, a/Rs, and l. We set the priors of B and b for both epochs, and Rp/Rs as uniform functions. We then set the priors of Tc for each epoch, u1, u2, and α for each band, Porb, a/Rs, and l as Gaussian priors. For the values and widths of the Gaussian priors, we referred to the values and widths of Gaussian priors from Johnson et al. (2015) for α and l and von Essen et al. (2014) for Porb; the others were obtained in the same manner as the spectral analysis. We set the logarithm of the likelihood ln Plike as
(3)
because we adopted Gaussian process for this fitting (Rasmussen & Williams 2005). Here, r is a series of residuals obtained by subtracting the model data from the observation data. In the MCMC process in PYMC (Salvatier, Wiecki & Fonnesbeck 2016), we ran 1000 steps, cut-off the first 500 steps as burn-in, and iterated this set 20 times. Figs A3 and A4 in Appendix  A plot the posterior distributions.

4 FITTING WITH NODAL PRECESSION MODEL

The angular momentum of WASP-33b’s planetary orbit |$|\overrightarrow{L_{p}}|$| (= 2πMpa2/Porb) is significantly smaller than the stellar rotational angular momentum of its host star |$|\overrightarrow{L_{s}}|$|⁠; |$|\overrightarrow{L_{p}}|/|\overrightarrow{L_{s}}|$| is ∼0.05 using the values of |$|\overrightarrow{L_{s}}|$| from Iorio (2011), those of Mp from Lehmann et al. (2015), a, and P from Collier Cameron et al. (2010). In this case, we can regard the stellar rotational axis as a stable vector and calculate the changes in b and λ
(4)
(5)
Here, θ, the nodal angle, can be expressed as follows:
(6)
where the slope of equation (6) is the precession speed from Barnes et al. (2013).
We then fitted the model in equations (4) and (5) with the values measured by the MCMC using PYMC. We considered ψ, θ(t = 2008), is, and J2 as free parameters and set their priors as uniform functions. Here, we set θ in 2008 as θ0, i.e. the initial value of θ in equation (6). For this fitting, we set the logarithm of the likelihood ln Plike as:
(7)
where Omes,i is the measured value of λ and b of each epoch, Cmod,i is the model value of λ and b, and σmes,i is the uncertainty of the measured λ and b, respectively. We consider that the values b in 2017 and 2018 were zero, as shown in Fig. A3. We ran 20 000 steps, cut-off the first 10 000 steps as burn-in, and iterated this set 20 times.
Table 1.

Details on the spectral data sets.

Date (UT)InstrumentNumber of spectraExposure time (s)ResolutionSNR at 5500 ÅReference
12 Nov 2008HJST/TS231390060 000140aCollier Cameron et al. (2010)
19 Oct 2011Subaru/HDS35600 (33 spectra), 480 (2 spectra)110 000160aWatanabe et al. (2020)
4 Oct 2014HJST/TS232190060 000280aJohnson et al. (2015)
28 Sep 2016TNG/HARPS-N40600115 000110bBorsa et al. (2021)
11 Dec 2016HJST/TS232190060 000250aThis work
12 Jan 2018TNG/HARPS-N23900115 000170bBorsa et al. (2021)
2 Jan 2019TNG/HARPS-N33600115 000120bBorsa et al. (2021)
27 Dec 2019OAC/MuSCAT12120065 00050aThis work
Date (UT)InstrumentNumber of spectraExposure time (s)ResolutionSNR at 5500 ÅReference
12 Nov 2008HJST/TS231390060 000140aCollier Cameron et al. (2010)
19 Oct 2011Subaru/HDS35600 (33 spectra), 480 (2 spectra)110 000160aWatanabe et al. (2020)
4 Oct 2014HJST/TS232190060 000280aJohnson et al. (2015)
28 Sep 2016TNG/HARPS-N40600115 000110bBorsa et al. (2021)
11 Dec 2016HJST/TS232190060 000250aThis work
12 Jan 2018TNG/HARPS-N23900115 000170bBorsa et al. (2021)
2 Jan 2019TNG/HARPS-N33600115 000120bBorsa et al. (2021)
27 Dec 2019OAC/MuSCAT12120065 00050aThis work
a

SNR per pixel.

b

SNR per extracted pixel.

Table 1.

Details on the spectral data sets.

Date (UT)InstrumentNumber of spectraExposure time (s)ResolutionSNR at 5500 ÅReference
12 Nov 2008HJST/TS231390060 000140aCollier Cameron et al. (2010)
19 Oct 2011Subaru/HDS35600 (33 spectra), 480 (2 spectra)110 000160aWatanabe et al. (2020)
4 Oct 2014HJST/TS232190060 000280aJohnson et al. (2015)
28 Sep 2016TNG/HARPS-N40600115 000110bBorsa et al. (2021)
11 Dec 2016HJST/TS232190060 000250aThis work
12 Jan 2018TNG/HARPS-N23900115 000170bBorsa et al. (2021)
2 Jan 2019TNG/HARPS-N33600115 000120bBorsa et al. (2021)
27 Dec 2019OAC/MuSCAT12120065 00050aThis work
Date (UT)InstrumentNumber of spectraExposure time (s)ResolutionSNR at 5500 ÅReference
12 Nov 2008HJST/TS231390060 000140aCollier Cameron et al. (2010)
19 Oct 2011Subaru/HDS35600 (33 spectra), 480 (2 spectra)110 000160aWatanabe et al. (2020)
4 Oct 2014HJST/TS232190060 000280aJohnson et al. (2015)
28 Sep 2016TNG/HARPS-N40600115 000110bBorsa et al. (2021)
11 Dec 2016HJST/TS232190060 000250aThis work
12 Jan 2018TNG/HARPS-N23900115 000170bBorsa et al. (2021)
2 Jan 2019TNG/HARPS-N33600115 000120bBorsa et al. (2021)
27 Dec 2019OAC/MuSCAT12120065 00050aThis work
a

SNR per pixel.

b

SNR per extracted pixel.

Table 2.

Measured parameters for WASP-33b.

Date (UT)λ (deg)bTc(BJDTDB)MethodInstrument
12 Nov 2008|$-111.30^{+0.76}_{-0.77}$||$0.2398^{+0.0062}_{-0.0058}$|2454782.92502 ± 0.00016Doppler tomographyHJST/TS23
19 Oct 2011−113.96 ± 0.300.1578 ± 0.0027|$2455853.96863^{+0.00014}_{-0.00015}$|Doppler tomographySubaru/HDS
4 Oct 2014−113.00 ± 0.370.0845 ± +0.00312456934.77139 ± 0.00015Doppler tomographyHJST/TS23
28 Sep 2016−111.39 ± 0.230.0413 ± 0.00192457660.59250 ± 0.00014Doppler tomographyTNG/HARPS-N
11 Dec 2016|$-111.32^{+0.49}_{-0.47}$|0.0432 ± 0.00392457733.78452 ± 0.00015Doppler tomographyHJST/TS23
5 Nov 2017|b| < 0.1322458063.14903 ± 0.00015PhotometryOAC/MuSCAT
11 Oct 2018|b| < 0.0672458403.49234 ± 0.00014PhotometryTCS/MuSCAT2
12 Jan 2018−111.46 ± 0.28|$0.0034^{+0.0024}_{-0.0023}$|2458131.46135 ± 0.00015Doppler tomographyTNG/HARPS-N
2 Nov 2019−111.64 ± 0.28|$-0.0272^{+0.0020}_{-0.0021}$|2458486.44296 ± 0.00016Doppler tomographyTNG/HARPS-N
27 Dec 2019|$-112.24^{+0.97}_{-1.02}$||$-0.0592^{+0.0066}_{-0.0065}$||$2458845.08379^{+0.00015}_{-0.00016}$|Doppler tomographyOAC/HIDES
Date (UT)λ (deg)bTc(BJDTDB)MethodInstrument
12 Nov 2008|$-111.30^{+0.76}_{-0.77}$||$0.2398^{+0.0062}_{-0.0058}$|2454782.92502 ± 0.00016Doppler tomographyHJST/TS23
19 Oct 2011−113.96 ± 0.300.1578 ± 0.0027|$2455853.96863^{+0.00014}_{-0.00015}$|Doppler tomographySubaru/HDS
4 Oct 2014−113.00 ± 0.370.0845 ± +0.00312456934.77139 ± 0.00015Doppler tomographyHJST/TS23
28 Sep 2016−111.39 ± 0.230.0413 ± 0.00192457660.59250 ± 0.00014Doppler tomographyTNG/HARPS-N
11 Dec 2016|$-111.32^{+0.49}_{-0.47}$|0.0432 ± 0.00392457733.78452 ± 0.00015Doppler tomographyHJST/TS23
5 Nov 2017|b| < 0.1322458063.14903 ± 0.00015PhotometryOAC/MuSCAT
11 Oct 2018|b| < 0.0672458403.49234 ± 0.00014PhotometryTCS/MuSCAT2
12 Jan 2018−111.46 ± 0.28|$0.0034^{+0.0024}_{-0.0023}$|2458131.46135 ± 0.00015Doppler tomographyTNG/HARPS-N
2 Nov 2019−111.64 ± 0.28|$-0.0272^{+0.0020}_{-0.0021}$|2458486.44296 ± 0.00016Doppler tomographyTNG/HARPS-N
27 Dec 2019|$-112.24^{+0.97}_{-1.02}$||$-0.0592^{+0.0066}_{-0.0065}$||$2458845.08379^{+0.00015}_{-0.00016}$|Doppler tomographyOAC/HIDES
Table 2.

Measured parameters for WASP-33b.

Date (UT)λ (deg)bTc(BJDTDB)MethodInstrument
12 Nov 2008|$-111.30^{+0.76}_{-0.77}$||$0.2398^{+0.0062}_{-0.0058}$|2454782.92502 ± 0.00016Doppler tomographyHJST/TS23
19 Oct 2011−113.96 ± 0.300.1578 ± 0.0027|$2455853.96863^{+0.00014}_{-0.00015}$|Doppler tomographySubaru/HDS
4 Oct 2014−113.00 ± 0.370.0845 ± +0.00312456934.77139 ± 0.00015Doppler tomographyHJST/TS23
28 Sep 2016−111.39 ± 0.230.0413 ± 0.00192457660.59250 ± 0.00014Doppler tomographyTNG/HARPS-N
11 Dec 2016|$-111.32^{+0.49}_{-0.47}$|0.0432 ± 0.00392457733.78452 ± 0.00015Doppler tomographyHJST/TS23
5 Nov 2017|b| < 0.1322458063.14903 ± 0.00015PhotometryOAC/MuSCAT
11 Oct 2018|b| < 0.0672458403.49234 ± 0.00014PhotometryTCS/MuSCAT2
12 Jan 2018−111.46 ± 0.28|$0.0034^{+0.0024}_{-0.0023}$|2458131.46135 ± 0.00015Doppler tomographyTNG/HARPS-N
2 Nov 2019−111.64 ± 0.28|$-0.0272^{+0.0020}_{-0.0021}$|2458486.44296 ± 0.00016Doppler tomographyTNG/HARPS-N
27 Dec 2019|$-112.24^{+0.97}_{-1.02}$||$-0.0592^{+0.0066}_{-0.0065}$||$2458845.08379^{+0.00015}_{-0.00016}$|Doppler tomographyOAC/HIDES
Date (UT)λ (deg)bTc(BJDTDB)MethodInstrument
12 Nov 2008|$-111.30^{+0.76}_{-0.77}$||$0.2398^{+0.0062}_{-0.0058}$|2454782.92502 ± 0.00016Doppler tomographyHJST/TS23
19 Oct 2011−113.96 ± 0.300.1578 ± 0.0027|$2455853.96863^{+0.00014}_{-0.00015}$|Doppler tomographySubaru/HDS
4 Oct 2014−113.00 ± 0.370.0845 ± +0.00312456934.77139 ± 0.00015Doppler tomographyHJST/TS23
28 Sep 2016−111.39 ± 0.230.0413 ± 0.00192457660.59250 ± 0.00014Doppler tomographyTNG/HARPS-N
11 Dec 2016|$-111.32^{+0.49}_{-0.47}$|0.0432 ± 0.00392457733.78452 ± 0.00015Doppler tomographyHJST/TS23
5 Nov 2017|b| < 0.1322458063.14903 ± 0.00015PhotometryOAC/MuSCAT
11 Oct 2018|b| < 0.0672458403.49234 ± 0.00014PhotometryTCS/MuSCAT2
12 Jan 2018−111.46 ± 0.28|$0.0034^{+0.0024}_{-0.0023}$|2458131.46135 ± 0.00015Doppler tomographyTNG/HARPS-N
2 Nov 2019−111.64 ± 0.28|$-0.0272^{+0.0020}_{-0.0021}$|2458486.44296 ± 0.00016Doppler tomographyTNG/HARPS-N
27 Dec 2019|$-112.24^{+0.97}_{-1.02}$||$-0.0592^{+0.0066}_{-0.0065}$||$2458845.08379^{+0.00015}_{-0.00016}$|Doppler tomographyOAC/HIDES
Table 3.

Calculated parameters for WASP-33b.

Dateψ (deg)J2is (deg)θ2008 (deg)
Results of this study|$108.19^{+0.95}_{-0.97}$||$(1.36^{+0.15}_{-0.12})\times 10^{-4}$||$58.3^{+4.6}_{-4.2}$||$73.6^{+1.7}_{-1.5}$|
Iorio (2016)|$99^{+5}_{-4}$| (in 2008), |$103^{+5}_{-4}$| (in 2014)|$(2.1^{+0.8}_{-0.5})\times 10^{-4}$||$142^{+10}_{-11}$|
Watanabe et al. (2020)(9.14 ± 0.51) × 10−5|$96^{+10}_{-14}$|
Borsa et al. (2021)113.99 ± 0.22(6.73 ± 0.22) × 10−590.11 ± 0.12
Dholakia, Luger & Dholakia (2022)|$108.3^{+19.0}_{-15.4}$||$69.8^{+4.0}_{-3.2}$|
Dateψ (deg)J2is (deg)θ2008 (deg)
Results of this study|$108.19^{+0.95}_{-0.97}$||$(1.36^{+0.15}_{-0.12})\times 10^{-4}$||$58.3^{+4.6}_{-4.2}$||$73.6^{+1.7}_{-1.5}$|
Iorio (2016)|$99^{+5}_{-4}$| (in 2008), |$103^{+5}_{-4}$| (in 2014)|$(2.1^{+0.8}_{-0.5})\times 10^{-4}$||$142^{+10}_{-11}$|
Watanabe et al. (2020)(9.14 ± 0.51) × 10−5|$96^{+10}_{-14}$|
Borsa et al. (2021)113.99 ± 0.22(6.73 ± 0.22) × 10−590.11 ± 0.12
Dholakia, Luger & Dholakia (2022)|$108.3^{+19.0}_{-15.4}$||$69.8^{+4.0}_{-3.2}$|
Table 3.

Calculated parameters for WASP-33b.

Dateψ (deg)J2is (deg)θ2008 (deg)
Results of this study|$108.19^{+0.95}_{-0.97}$||$(1.36^{+0.15}_{-0.12})\times 10^{-4}$||$58.3^{+4.6}_{-4.2}$||$73.6^{+1.7}_{-1.5}$|
Iorio (2016)|$99^{+5}_{-4}$| (in 2008), |$103^{+5}_{-4}$| (in 2014)|$(2.1^{+0.8}_{-0.5})\times 10^{-4}$||$142^{+10}_{-11}$|
Watanabe et al. (2020)(9.14 ± 0.51) × 10−5|$96^{+10}_{-14}$|
Borsa et al. (2021)113.99 ± 0.22(6.73 ± 0.22) × 10−590.11 ± 0.12
Dholakia, Luger & Dholakia (2022)|$108.3^{+19.0}_{-15.4}$||$69.8^{+4.0}_{-3.2}$|
Dateψ (deg)J2is (deg)θ2008 (deg)
Results of this study|$108.19^{+0.95}_{-0.97}$||$(1.36^{+0.15}_{-0.12})\times 10^{-4}$||$58.3^{+4.6}_{-4.2}$||$73.6^{+1.7}_{-1.5}$|
Iorio (2016)|$99^{+5}_{-4}$| (in 2008), |$103^{+5}_{-4}$| (in 2014)|$(2.1^{+0.8}_{-0.5})\times 10^{-4}$||$142^{+10}_{-11}$|
Watanabe et al. (2020)(9.14 ± 0.51) × 10−5|$96^{+10}_{-14}$|
Borsa et al. (2021)113.99 ± 0.22(6.73 ± 0.22) × 10−590.11 ± 0.12
Dholakia, Luger & Dholakia (2022)|$108.3^{+19.0}_{-15.4}$||$69.8^{+4.0}_{-3.2}$|

5 RESULTS

Fig. 2 shows the line profile residuals and the best-fitted filtered models. Table 2 lists the best values of λ and b. Our results for λ and b in 2014 were in excellent agreement with the values of Johnson et al. (2015) within 1σ, whereas those in 2008 were marginally consistent with Johnson et al. (2015) within 2σ. Moreover, our results for λ in 2016 and 2018 from HARPS-N were consistent with those of Borsa et al. (2021), whereas our λ in 2019 from HARPS-N differs from that in Borsa et al. (2021) by ∼2σ.

Figs 3 and 4 show the best-fitting light-curve models. Table 2 lists the ranges of b from MuSCAT and MuSCAT2. The posteriors of both b in Fig. A3 exhibit a truncated normal distribution with a minimum value of ∼ 0. Hence, we set 1σ as a 68 per cent confidence interval based on the minimum value of 0 because they were identical, but with opposite signs of impact parameters, which yielded the same transit light curves.

Moreover, Fig. 5 presents posteriors from MCMC with the nodal precession model. Table 3 lists the values of ψ, θ(t = 2008), is, and J2. Here, we note that Iorio (2016) may have considered ψ as a variable value. Fig. 6 shows the change in λ and b for WASP-33b.

6 DISCUSSION

We have inspected the nodal precession of WASP-33b with more observations than that used in previous studies (Johnson et al. 2015; Watanabe et al. 2020; Borsa et al. 2021). Our study is the first to verify the nodal precession based on both of Doppler tomographic observation and transit photometry. The errors from the transit photometric observations were large. However, with the change in the impact parameter (see the left bottom part in Fig. 6), the results from the transit photometry were consistent with the predicted values from the decreasing trend for the Doppler tomographic observations. This indicates that transit photometry can simultaneously contribute to the measurements of the nodal precession using Doppler tomographic data. In Fig. 6, although the impact parameter of WASP-33b appeared to change linearly, the change in λ may not be along the model of the nodal precession. Thus, we should observe WASP-33b via Doppler tomography to clarify whether its λ increases based on the model or decreases from 2021.

6.1 Comparison of stellar spin inclination and quadrupole moment with previous studies

Our is value disagrees with those in previous studies (Johnson et al. 2015; Iorio 2016; Watanabe et al. 2020; Borsa et al. 2021), as derived from the nodal precession by ∼3σ or more. However, our value is similar to that of Dholakia et al. (2022), despite values that are ∼2σ different. Notably, Dholakia et al. (2022) derived WASP-33b’s is from light curve TESS photometric data considering its oblateness and gravity darkening.

The derived stellar quadrupole moment of WASP-33b was |$J_2 = (1.36^{+0.15}_{-0.12}) \times 10^{-4}$|⁠. This value agrees with that of Iorio (2016) within 1.5σ, which is larger than those of the other previous studies (add> 3σ) and smaller than the theoretical value (J2 = 3.8 × 10−4) calculated by Iorio (2011).

One of the possible reasons for disagreements with the values of is and J2 could be the difference in the nodal precession model. Iorio (2016) and Watanabe et al. (2020) used time variation models for other orbital parameters: the ascending node Ω and the orbital inclination to the apparent equatorial plane I, illustrated in fig. 4 of Watanabe et al. (2020), and calculated from λ and b. Then, Borsa et al. (2021) estimated the change in the inclination angle |$i_p (=\arccos {(bR_s/a)})$| as a linear function. In this study, we directly used the accurate time variation models of λ and b, such that our derived values for is and J2 would be more accurate than those in previous studies. However, we should clarify the cause of the short-term variation of λ to create more detailed nodal precession model.

In contrast, there is a probability that adding data sets also causes the disagreements. In this study, we found a short-term variation in λ, which may have been decreasing since 2016, although the reason remains unclear. Therefore, we have to obtain more data sets of WASP-33b’s transits to more accurately determine their values.

6.2 Orbital evolution of WASP-33b

We found the real spin–orbit obliquity of WASP-33b |$\psi =108.19^{+0.95}_{-0.97}$| deg. This uncertainty is larger than that reported by Borsa et al. (2021) owing to the difference in the precession model, which shows that dλ/dt of WASP-33b is always positive, while that of Borsa et al. (2021) allows dλ/dt = 0 at a certain time. The derived value indicates the possibility that WASP-33b has experienced the planet–planet scattering (Chatterjee et al. 2008) or Kozai migration (Fabrycky & Tremaine 2007), which are the mechanisms that cause the misaligned orbit. The existence of the WASP-33b’s companion star is necessary for clarifying and distinguishing these two evolutionary models. Nevertheless, Doppler tomographic observations and transit photometry cannot detect the companion stars.

Ngo et al. (2016) found a companion candidate for WASP-33, estimated as a dwarf star or brown dwarf and located at 238 au (Porb ∼ 3300 yr) from the host star based on direct imaging. However, this candidate has not yet been confirmed owing to negligible proper motion. Thus, additional direct imaging observations are required to confirm whether the companion moves in the same proper motion as the host star or not. If the stellar companion candidate, with 0.1 M, revolves in a circular orbit (e = 0) and WASP-33b was formed near the snow line, ∼5 au (Porb ∼ 12 yr), the planet experienced a Kozai oscillation with the Kozai cycle period, PKozai, of 13 Myr based on the following expression:
(8)
where Ms, Mc, Pc, Pb, and ec are the stellar mass, companion’s mass, companion’s period, planetary mass, and companion’s eccentricity, respectively (Wu, Murray & Ramsahai 2007). Although we have to consider that the distance from the host star to the companion star is the sky-projected distance and the eccentricity of the companion star may exist, the Kozai cycle period could be shorter than the age of WASP-33 (∼ 100 Myr). Therefore, the stellar companion may have caused the Kozai mechanism for WASP-33b.

6.3 Nodal precession speed

We calculated the nodal precession speed |$\dot{\theta }=0.507^{+0.025}_{-0.022}$| deg yr−1 and its period |$P_{\mathrm{pre}}=709^{+33}_{-34}$| yr, where the precession was faster than that of Iorio (2016) (⁠|$\dot{\theta }=0.37^{+0.04}_{-0.03}$|⁠). We then found that WASP-33b transits in front of the host star for only ∼ 20 per cent of the entire nodal precession period, which indicates that it is rare to discover WASP-33b as a transiting planet. This implies that WASP-33b began transiting in 1977 ± 2 and will stop transiting in 2055 ± 2.

6.4 Nodal precession observations of other hot Jupiters around hot stars

Doppler tomography has confirmed 17 hot Jupiters around hot stars to date. The real spin–orbit obliquities of seven of these planets (WASP-33b, Kepler-13Ab, KELT-9b, KELT-17b, MASCARA-1b, MASCARA-4b, and WASP-189b) have been identified thus far, which are nearly vertical (60 deg < ψ < 120 deg). The remaining hot Jupiters reveal that only their projected spin–orbit obliquities were obtained. Three hot Jupiters around hot stars with λ ∼ 90 deg, i.e. KELT-26b (⁠|$\lambda =91.3^{+6.5}_{-6.3}$| deg; Rodríguez Martínez et al. 2020), HAT-P-70b (⁠|$\lambda =113.1^{+5.1}_{-3.4}$| deg; Zhou et al. 2019), and TOI-1518b (⁠|$\lambda =-119.66^{+0.98}_{-0.93}$| deg; Cabot et al. 2021), are valuable for detecting their observable nodal precessions to derive their ψ. Even if λ is near 0 deg or 180 deg, we can observe the change in the transit trajectory and measure ψ near 90 deg when the star rotates nearly pole-on for the line of sight. However, when λ is near 0 or 180 deg and the star rotation axis is almost perpendicular to the line of sight, the transit trajectory barely moves because ψ should also be near 0 deg or 180 deg. In this case, we can estimate ψ as λ. Therefore, we should regularly observe the nodal precessions of these planets around their hot stars to measure their ψ using Doppler tomography and transit photometry.

Kepler-13Ab, orbiting around an A-type star, is the another hot Jupiter whose nodal precession has been detected by only transit photometries (Barnes, Linscott & Shporer 2011; Szabó et al. 2011). Previous studies have measured ψ of this hot Jupiter via gravity-darkened transit photometry, but the values were different between these two results (ψ = 60 ± 2 deg in Masuda 2015 and ψ = 29 ± 1 deg in Herman et al. 2018). Although Kepler-13Ab is likely to have evolved with Kozai migration owing to its companion star, Kepler-13B (Santerne et al. 2012), the value of its ψ should be verified by adding Doppler tomographic observations and transit photometries to obtain an accurate histogram of the derived ψ. As Johnson et al. (2014) measured the λ of Kepler-13Ab in 2014 via Doppler tomography, an additional transit spectroscopic observation enabled us to detect the change in λ and then to derive its ψ independently from the gravity-darkened transit photometry.

Ahlers et al. (2020a, b), Lendl et al. (2020), and Hooton et al. (2022) measured the ψ of MASCARA-4b (⁠|$\psi =104^{+7}_{-13}$| deg), WASP-189b (ψ = 85.4 ± 4.3 deg), KELT-9b (⁠|$\psi =87^{+10}_{-11}$| deg), and MASCARA-1b (⁠|$\psi =72.1^{+2.5}_{-2.4}$| deg), respectively, using gravity-darkened transit photometry. Zhou et al. (2016) derived the ψ of KLET-17b (ψ = 116 ± 4 deg) by the technique using the differential rotation of a host star.

These planets are hot Jupiters around A-type stars: their nodal precessions are yet to be detected. Observing the nodal precessions of these hot Jupiters is important for verifying the values of their ψ. This observation can also contribute to the investigations on the internal structure of hot stars by deriving the values of their J2.

Although there are 17 hot Jupiters around hot stars whose projected spin–orbit obliquities have been measured, this number remains too small to statistically determine the orbital evolution tendency. Albrecht et al. (2021) found that planets around solar-like stars with large projected spin–orbit obliquities are likely to revolve on polar orbits (ψ ∼ 90 deg); however, the tendency remains unclear for hot Jupiters around hot stars because the real spin–orbit obliquities of only six of them have been revealed. If normal planet–disc interaction is the main migration, the distribution should gather at ψ = 0 deg (Lai, Foucart & Lin 2011). If the orbital evolution by planet–planet scattering is the majority, the orbits are likely to incline at approximately ψ = 60 deg (Nagasawa & Ida 2011). When Kozai migration is the primary evolution, ψ has a wide range from 10 deg to 140 deg (Petrovich 2015). As Barclay, Pepper & Quintana (2018) predicted that the TESS mission could find 500 hot Jupiters around A-type stars from the 2-yr observation, increasing the number of hot Jupiters around hot stars will become possible by validating TESS planet candidates. Further observations could lead to uncovering the detailed ψ distribution around hot stars.

7 CONCLUSION

We analysed the nodal precession of WASP-33b via Doppler tomography and transit photometry using various high-dispersion spectrographs including Subaru/HDS, HJST/TS23, OAC/HIDES, and TNG/HARPS-N, and two multicolour simultaneous cameras, OAC/MuSCAT and TCS/MuSCAT2. Based on the observed change in the projected spin–orbit obliquity λ and the impact parameter b, we modelled the nodal precession of WASP-33b and derived the real spin–orbit obliquity of WASP-33b as |$\psi = 108.19^{+0.95}_{-0.97}$| deg. Compared with the results of previous studies, the results of the near-polar orbit did not change. However, our value for ψ differed from that of Borsa et al. (2021) by >3σ. This discrepancy may be caused by whether or not psi is a constant value. We assumed psi as a constant while Borsa et al. (2021) did not; Borsa et al. (2021) adopted psi from its value at the 2011 epoch. We also simultaneously derived the stellar spin inclination and the stellar gravitational quadrupole moment of WASP-33 as |$i_\mathit{ s} = 58.3^{+4.6}_{-4.2}$| deg, and |$J_2=(1.36^{+0.15}_{-0.12}) \times 10^{-4}$|⁠, respectively. These results differed by >3σ from those of previous studies on nodal precession, except for the J2 value reported in Iorio (2016) (∼2.1σ).

A likely reason for these discrepancies is the different nodal precession models. For the first time, we applied the accurate time variation models for λ and b to fit the nodal precession. Therefore, our derived values of ψ, is, and J2 are the most accurate to date. Moreover, additional data sets may have updated these parameter values, thus causing these differences. Thus, acquiring more data sets will allow for the derivation of more accurate values of ψ, is, and J2.

We calculated the nodal precession speed |$\dot{\theta }=0.507^{+0.025}_{-0.022}$| deg yr−1 and its period |$P_{\mathrm{pre}}=709^{+33}_{-34}$| yr, which revealed that WASP-33b transits in front of the host star for only ∼20 per cent of the entire nodal precession period. Based on this result, we speculate that WASP-33b began transiting in 1977 ± 2; we forecast that it will finish transiting in 2055 ± 2.

The TESS survey should help us increase our ability to count number of hot Jupiters around hot stars in the future. Applying the proposed methodology to newly discovered hot Jupiters around hot stars is important not only for characterizing each planetary system, but also for discriminating the migration mechanisms of such planets and investigating the internal structure of hot stars.

ACKNOWLEDGEMENTS

This paper is based on data collected at the Subaru Telescope, which is located atop Maunakea and operated by the National Astronomical Observatory of Japan (NAOJ). We wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. The paper also includes data taken at The McDonald Observatory of The University of Texas at Austin and taken at The Okayama Astrophysical Observatory. This article is based on observations made with the MuSCAT2 instrument, developed by Astrobiology Center (ABC), at Telescopio Carlos Sánchez operated on the island of Tenerife by the Instituto de Astrofísica de Canarias (IAC) in the Spanish Observatorio del Teide. pyraf is a product of the Space Telescope Science Institute, which is operated by Association of Universities for Research in Astronomy (AURA) for National Aeronautics and Space Administration (NASA). This work has made use of the Vienna Atomic Line Database (VALD) data base, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. We are grateful to editage (https://www.editage.jp) for English editing. We acknowledge the Global Architecture of Planetary Systems (GAPS) Consortium (Covino et al. 2013) for providing the mean line profiles of their HARPS-N transits, and thank F. Borsa for providing spectral data sets of HARPS-N. This work is partly supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers JP21K20376, JP18H05439, JP20J21872, JP17H04574,JP20K14518, Japan Science and Technology Agency (JST) CREST Grant Number JPMJCR1761, Astrobiology Center SATELLITE Research project AB022006, and the Astrobiology Center of National Institutes of Natural Sciences (NINS) (Grant Number AB031010). N. B-C. acknowledges funding from the European Research Council under the European Union’s Horizon 2020 research and innovation program under grant agreement number 694513. E. E-B. acknowledges financial support from the European Union and the State Agency of Investigation of the Spanish Ministry of Science and Innovation (MICINN) under the grant PRE2020-093107 of the Pre-Doc Program for the Training of Doctors (FPI-SO) through FSE funds.

DATA AVAILABILITY

The raw data from OAC/HIDES, OAC/MuSCAT, and TCS/MuSCAT2 will be shared on reasonable request to the corresponding author. The reduced data of Subaru/HDS were provided by N. Narita by permission, and will be shared on request to the corresponding author with permission of N. Narita. The line profile data from HJST/TS23 were provided by M. C. Johnson by permission, and will be shared on request to the corresponding author with permission of M. C. Johnson. The line profile data of TNG/HARPS-N were provided by F. Borsa by permission, and will be shared on request to the corresponding author with permission of F. Borsa.

REFERENCES

Abt
H. A.
,
Levato
H.
,
Grosso
M.
,
2002
,
ApJ
,
573
,
359

Ahlers
J. P.
et al. ,
2020a
,
AJ
,
160
,
4

Ahlers
J. P.
et al. ,
2020b
,
ApJ
,
888
,
63

Albrecht
S.
et al. ,
2012
,
ApJ
,
757
,
18

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

Barclay
T.
,
Pepper
J.
,
Quintana
E. V.
,
2018
,
ApJS
,
239
,
2

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

Barnes
J. W.
,
van Eyken
J. C.
,
Jackson
B. K.
,
Ciardi
D. R.
,
Fortney
J. J.
,
2013
,
ApJ
,
774
,
53

Borsa
F.
et al. ,
2021
,
A&A
,
653
,
A104

Cabot
S. H. C.
et al. ,
2021
,
AJ
,
162
,
218

Chatterjee
S.
,
Ford
E. B.
,
Matsumura
S.
,
Rasio
F. A.
,
2008
,
ApJ
,
686
,
580

Collier Cameron
A.
et al. ,
2007
,
MNRAS
,
375
,
951

Collier Cameron
A.
et al. ,
2010
,
MNRAS
,
407
,
507

Cosentino
R.
et al. ,
2012
, in
McLean
I. S.
,
Ramsay
S. K.
,
Takami
H.
, eds,
Proc. SPIE Conf. Ser. Vol. 8446, Ground-Based and Airborne Instrumentation for Astronomy IV
.
SPIE
,
Bellingham
, p.
84461V

Covino
 
E.
et al. ,
2013
,
A&A
,
554
,
A28

Dholakia
S.
,
Luger
R.
,
Dholakia
S.
,
2022
,
ApJ
,
925
,
185

Donati
J.-F.
,
Semel
M.
,
Carter
B. D.
,
Rees
D. E.
,
Collier Cameron
A.
,
1997
,
MNRAS
,
291
,
658

Fabrycky
D.
,
Tremaine
S.
,
2007
,
ApJ
,
669
,
1298

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

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

Foreman-Mackey
D.
,
Luger
R.
,
Czekala
I.
,
Agol
E.
,
Price-Whelan
A.
,
Brandt
T. D.
,
Barclay
T.
,
Bouma
L.
,
2020
,
exoplanet-dev/exoplanet v0.4.0
.

Fukui
A.
et al. ,
2011
,
PASJ
,
63
,
287

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

Herrero
E.
,
Morales
J. C.
,
Ribas
I.
,
Naves
R.
,
2011
,
A&A
,
526
,
L10

Hooton
M. J.
et al. ,
2022
,
A&A
,
658
,
A75

Husser
T.-O.
,
Wende-von Berg
S.
,
Dreizler
S.
,
Homeier
D.
,
Reiners
A.
,
Barman
T.
,
Hauschildt
P. H.
,
2013
,
A&A
,
553
,
A6

Iorio
L.
,
2011
,
Ap&SS
,
331
,
485

Iorio
L.
,
2016
,
MNRAS
,
455
,
207

Izumiura
H.
,
1999
, in
Chen
P. S.
, ed.,
Proc. 4th East Asian Meeting on Astronomy, Observational Astrophysics in Asia and its Future
. p.
77

Johnson
M. C.
,
Cochran
W. D.
,
Albrecht
S.
,
Dodson-Robinson
S. E.
,
Winn
J. N.
,
Gullikson
K.
,
2014
,
ApJ
,
790
,
30

Johnson
M. C.
,
Cochran
W. D.
,
Collier Cameron
A.
,
Bayliss
D.
,
2015
,
ApJ
,
810
,
L23

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

Kochukhov
O.
,
Makaganiuk
V.
,
Piskunov
N.
,
2010
,
A&A
,
524
,
A5

Kovács
G.
et al. ,
2013
,
A&A
,
553
,
A44

Kupka
F. G.
,
Ryabchikova
T. A.
,
Piskunov
N. E.
,
Stempels
H. C.
,
Weiss
W. W.
,
2000
,
Balt. Astron.
,
9
,
590

Lai
D.
,
Foucart
F.
,
Lin
D. N. C.
,
2011
,
MNRAS
,
412
,
2790

Lehmann
H.
,
Guenther
E.
,
Sebastian
D.
,
Döllinger
M.
,
Hartmann
M.
,
Mkrtichian
D. E.
,
2015
,
A&A
,
578
,
L4

Lendl
M.
et al. ,
2020
,
A&A
,
643
,
A94

Masuda
K.
,
2015
,
ApJ
,
805
,
28

Nagasawa
M.
,
Ida
S.
,
2011
,
ApJ
,
742
,
72

Narita
N.
et al. ,
2015
,
J. Astron. Telesc. Instrum. Syst.
,
1
,
045001

Narita
N.
et al. ,
2019
,
J. Astron. Telesc. Instrum. Syst.
,
5
,
015001

Ngo
H.
et al. ,
2016
,
ApJ
,
827
,
8

Noguchi
K.
et al. ,
2002
,
PASJ
,
54
,
855

Parviainen
H.
,
Aigrain
S.
,
2015
,
MNRAS
,
453
,
3821

Pepper
J.
et al. ,
2007
,
PASP
,
119
,
923

Petrovich
C.
,
2015
,
ApJ
,
799
,
27

Rasmussen
C. E.
,
Williams
C. K. I.
,
2005
,
Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)
.
The MIT Press
,
Cambridge, MA

Ricker
G. R.
et al. ,
2015
,
J. Astron. Telesc. Instrum. Syst.
,
1
,
014003

Rodríguez Martínez
R.
et al. ,
2020
,
AJ
,
160
,
111

Salvatier
J.
,
Wiecki
T. V.
,
Fonnesbeck
C.
,
2016
,
Peer J Comput. Sci.
,
2
,
e55

Santerne
A.
et al. ,
2012
,
A&A
,
544
,
L12

Science Software Branch at STScI
,
2012
,
Astrophysics Source Code Library
,
record ascl:1207.011

Szabó
G. M.
et al. ,
2011
,
ApJ
,
736
,
L4

Szabó
G. M.
,
Pál
A.
,
Derekas
A.
,
Simon
A. E.
,
Szalai
T.
,
Kiss
L. L.
,
2012
,
MNRAS
,
421
,
L122

Tull
R. G.
,
MacQueen
P. J.
,
Sneden
C.
,
Lambert
D. L.
,
1995
,
PASP
,
107
,
251

von Essen
C.
et al. ,
2014
,
A&A
,
561
,
A48

von Essen
C.
,
Mallonn
M.
,
Welbanks
L.
,
Madhusudhan
N.
,
Pinhas
A.
,
Bouy
H.
,
Weis Hansen
P.
,
2019
,
A&A
,
622
,
A71

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

Wu
Y.
,
Murray
N. W.
,
Ramsahai
J. M.
,
2007
,
ApJ
,
670
,
820

Zhou
G.
et al. ,
2016
,
AJ
,
152
,
136

Zhou
G.
et al. ,
2019
,
AJ
,
158
,
141

APPENDIX A: MCMC RESULTS OF DOPPLER TOMOGRAPHIC AND PHOTOMETRIC MEASUREMENTS

Here, we display the corner plots after using MCMC in Section 2.3 (Figs A1 and A2) and in Section 3.2 (Figs A3 and A4).

Corner plots for the free parameters after using MCMC in Section 2.3. Black circles indicate 68 per cent, 95 per cent, and 99.7 per cent confidence from the inside. In each posterior distribution of each parameter, the vertical dotted lines show the best-fitting value (middle) and 1σ confidence (both ends). We created these plots using corner.py (Foreman-Mackey 2016).
Figure A1.

Corner plots for the free parameters after using MCMC in Section 2.3. Black circles indicate 68 per cent, 95 per cent, and 99.7 per cent confidence from the inside. In each posterior distribution of each parameter, the vertical dotted lines show the best-fitting value (middle) and 1σ confidence (both ends). We created these plots using corner.py (Foreman-Mackey 2016).

Continuance of Fig. A1.
Figure A2.

Continuance of Fig. A1.

Plot distributions for impact parameter of each epoch, radial-ratio, period and transit mid-time of each epoch, limb darkening coefficients of each band, and base line of each light curve. These are the same corner plots as in Fig. A1, but for the transit photometry by MuSCAT1 and MuSCAT2.
Figure A3.

Plot distributions for impact parameter of each epoch, radial-ratio, period and transit mid-time of each epoch, limb darkening coefficients of each band, and base line of each light curve. These are the same corner plots as in Fig. A1, but for the transit photometry by MuSCAT1 and MuSCAT2.

Continuance of Fig. A3. Plot distributions for the impact parameter of each epoch, σ of each band, l, period, and transit mid-time of each epoch.
Figure A4.

Continuance of Fig. A3. Plot distributions for the impact parameter of each epoch, σ of each band, l, period, and transit mid-time of each epoch.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)