Rossiter–McLaughlin models and their effect on estimates of stellar rotation, illustrated using six WASP systems^{†}^{†}thanks: based on observations (under proposal 090.C0540) made using the HARPS high resolution échelle spectrograph mounted on the ESO 3.6 m at the ESO La Silla observatory, and completed by photometry obtained the Swiss 1.2m Euler Telescope, also at La Silla.
Abstract
We present new measurements of the projected spin–orbit angle for six WASP hot Jupiters, four of which are new to the literature (WASP61, 62, 76, and 78), and two of which are new analyses of previously measured systems using new data (WASP71, and 79). We use three different models based on two different techniques: radial velocity measurements of the Rossiter–McLaughlin effect, and Doppler tomography. Our comparison of the different models reveals that they produce projected stellar rotation velocities () measurements often in disagreement with each other and with estimates obtained from spectral line broadening. The Boué model for the Rossiter–McLaughlin effect consistently underestimates the value of compared to the Hirano model. Although differed, the effect on was small for our sample, with all three methods producing values in agreement with each other. Using Doppler tomography, we find that WASP61 b (), WASP71 b (), and WASP78 b () are aligned. WASP62 b () is found to be slightly misaligned, while WASP79 b () is confirmed to be strongly misaligned and has a retrograde orbit. We explore a range of possibilities for the orbit of WASP76 b, finding that the orbit is likely to be strongly misaligned in the positive direction.
keywords:
techniques: photometric – techniques: radial velocities – techniques: spectroscopic – planetary systems – stars: rotation1 Introduction
All eight planets of the Solar system orbit in approximately the same plane, the ecliptic, which is inclined to the solar equatorial plane by only (Beck & Giles, 2005). The orbital axes for the Solar system planets therefore exhibit near spinorbit alignment with the Sun’s rotation axis (the origin of the slight divergence from true alignment is unknown). There is no guarantee, however, that this holds true for extrasolar planets, as it is known that from binary stars that spinorbit angles can take a wide variety of values (e.g. Hube & Couch, 1982; Hale, 1994; Albrecht et al., 2007, 2009; Jensen & Akeson, 2014; Albrecht et al., 2014). Compared to such systems, measurement of the alignment angle (‘obliquity’) for an extrasolar planet is more difficult owing to the greater radius and luminosity ratios. This is compounded by the face that the host star of a closein exoplanet generally rotates more slowly than does the primary star in a stellar binary of the same orbital period. We are also generally limited to measuring the alignment angle as projected on to the plane of the sky, generally referred to as . Measurement of the true obliquity () requires knowledge of the inclination of the stellar rotation axis to the line of sight, , an angle that is currently very difficult to measure directly. It is possible to infer a value for using knowledge of the projected stellar rotation speed, , the stellar radius, , and the stellar rotation period, (e.g. Lund et al., 2014) , but the last of these can in turn be tricky to determine (e.g. Lendl et al., 2014).
HD 209458 b was the first extrasolar planet for which was measured (Queloz et al., 2000a). Since that work the number of systems for which the projected spinorbit alignment angle has been measured or inferred has been increasing at a steady rate, and is now closing in on 100. By far the majority of these are transiting, hot Jupiter extrextrasolarasolar planets, and for most of these the value of has been modelled using the Holt–Rossiter–McLaughlin (RM) effect (Holt, 1983; Schlesinger, 1910, 1916; Rossiter, 1924; McLaughlin, 1924), the spectroscopic signature that is produced during a transit by the occultation of the red and blueshifted stellar hemispheres. Other, complementary methods such as Doppler tomography (DT, Collier Cameron et al. 2010a), consideration of the gravity darkening effect (e.g. Barnes, Linscott & Shporer, 2011), modelling of photometric starspot signatures (e.g. Nutzman, Fabrycky & Fortney, 2011; SanchisOjeda et al., 2011; TregloanReed et al., 2015), measurement of the chromospheric RM effect in the Ca II H & K lines (Czesla et al., 2012), and analysis of photometric variability distributions (Mazeh et al., 2015) have also contributed to the tally.
While the vast majority of the measurements have been made via the RM effect, the models that have been used to model this effect have changed over time, becoming more complex and incorporating more detailed physics. The first models in widespread use were those of Ohta, Taruya & Suto (2005, 2009) and Giménez (2006), but these were superseded by the more detailed models of Hirano et al. (2011) and Boué et al. (2013), which take different approaches to the problem. A recent addition to the stable of RM models is that of Baluev & Shaidulin (2015). This assortment of models means, combined with the variety of instruments with which RV measurements are made, might be introducing biases into the parameters that we measure, particularly and . These have yet to be fully explored.
In this work, we present analysis of the spinorbit alignment in six hot Jupiter systems, found by the WASP consortium (Pollacco et al., 2006), with the aim of shedding new light on the problems discussed above. We observed WASP61 (Hellier et al., 2012); WASP62 (Hellier et al., 2012); WASP71 (Smith et al., 2013) WASP76 (West et al., 2016); WASP78 (Smalley et al., 2012), and WASP79 (Smalley et al., 2012).
These six systems were observed with HARPS under programme ID 090.C0540 (PI Triaud). Our earlier RM observation campaigns selected systems across a wide range of parameters and aimed to increase the number of spinorbit measurements with as few preconceptions as possible. Here, we instead selected six particular objects. At the time, Schlaufman (2010), Winn et al. (2010) and Triaud (2011) had noticed intriguing relations between some stellar parameters and the projected spinorbit angle. Our selection of planets, orbiting stars with around 6250 K, was meant to verify these.
2 Methods
As in Brown et al. (2012a, b), we analyse the complete set of available data for each system: WASP photometry; followup photometric transit data from previous studies; followup spectroscopic data from previous studies; newly acquired photometric transit data, and newly acquired intransit spectroscopic measurements of the RM effect using HARPS. New data are described in the appropriates subsections of Section 3, and are available both in the appendix and online as supplementary information. Our modelling has been extensively described in previous papers from the SuperWASP collaboration (Collier Cameron et al., 2007; Pollacco et al., 2008; Brown et al., 2012b, e.g.), but we summarize the process here for new readers.
Our analysis is carried out using a Markov Chain Monte Carlo (MCMC) algorithm using the Metropolis–Hastings decision maker (Metropolis et al., 1953; Hastings, 1970). Our jump parameters are listed in Table 1, and have been formulated to minimize correlations and maximize mutual orthogonality between parameters. We use , , , and to impose uniform priors on and and avoid bias towards higher values (Ford, 2006; Anderson et al., 2011). Several of our parameters (namely impact parameter, , [Fe/H], and the Boué model parameters when appropriate) are controlled by Gaussian priors by default (see Table 2). Others may be controlled by a prior if desired, see Section 2.2.
At each MCMC step, we calculate models of the photometric transit (following Mandel & Algol 2002), the Keplerian RV curve, and the RM effect (see following section). Photometric data is linearly decorrelated to remove systematic trends. Limbdarkening is accounted for by a fourcomponent, nonlinear model, with wavelength appropriate coefficients derived at each MCMC step by interpolation through the tables of Claret (2000, 2004). The RV Keplerian curve, and thus the orbital elements, is primarily constrained by the existing spectroscopic data, as our new, intransit spectroscopy covers only a small portion of the orbital phase. Quality of fit for these models is determined by calculating .
Other parameters are derived at each MCMC step using standard methods. Stellar mass, for example, is calculated using the calibration of Torres et al. (2010), with updated parameters from Southworth (2011). Stellar radius is calculated from (derived directly from the transit model) and the orbital period, via Kepler’s third law.
We use a burnin phase with a minimum of 500 steps, judging the chain to be converged (and thus burnin complete) when for that step is greater than the median of all previous values from the burnin chain (Knutson et al., 2008). This is followed by a phase of 100 accepted steps, which are used to rescale the error bars on the primary jump parameters, and a production run of accepted steps. Five separate chains are run, and the results concatenated to produce the final chain of length steps. The reported parameters are the median values from this final chain, with the uncertainties taken to be the values that enclose percent of the distribution.
We test for convergence of our chains using the statistics of (Geweke 1992, to check interchain convergence) and (Gelman & Rubin 1992, to check intrachain convergence). We also carry out additional visual checks using trace plots, autocorrelation plots, and probability distribution plots (both one and twodimensional). If an individual chain is found to be unconverged then we run a replacement chain, recalculating the reported parameters and convergence statistics. This process is repeated as necessary until the convergence tests indicate a fullyconverged final chain.
We have used the UTC time standard and Barycentric Julian Dates in our analysis. Our results are based on the equatorial solar and jovian radii, and masses, taken from Allen’s Astrophysical Quantities.
Parameter  Units  Symbol  Prior? 

Epoch  No  
Orbital period  days  No  
Transit width  days  No  
Transit depth  –  No  
Impact parameter  Stellar radii  Yes  
Effective temperature  K  Yes  
‘Metallicity’  dex  Yes  
RV semiamplitude  km s  No  
–  indirectly; Yes/No  
–  indirectly; Yes/No  
Longterm RV trend  –  Yes/No  
(km s)  indirectly; Yes/No  
(km s)  indirectly; Yes/No  
Barycentric RV for CCFs  km s  Yes/No  
FWHM for CCFs  km s  FWHM  No 
Barycentric RV for RM data  km s  Yes  
Boué model Gaussian width  km s  Yes 
System  Parameter  

WASP61  
WASP62  
WASP71  
WASP76  
WASP78  
WASP79 
2.1 Modelling spinorbit alignment
Our first model for the RM effect is that of Hirano et al. (2011). This has become the de facto standard thanks to its rigourous approach to the fitting procedure, which crosscorrelates an intransit spectrum with a template, and maximizes the crosscorrelation function (CCF). This method requires prior knowledge of several broadening coefficients, specifically the macroturbulence, , and the Lorentzian () and Gaussian () spectral line dispersions. For this work we assumed km s in line with Hirano et al., and also assumed that the coefficient of differential rotation, ^{*}^{*}*Whilst several of the systems under consideration are rapidly rotating, without knowledge of the inclination of their stellar rotation axes it is difficult to place a value on .. is calculated individually for each RV data set, and depends on the instrument used to collect the data as it is a function of the spectral resolution.
Boué et al. (2013) pointed out that the Hirano et al. model is poorly optimized for instruments which use a CCF based approach to their data reduction. For iodine cell spectrographs (e.g. HIRES at the Keck telescope), the Hirano et al. (2011) model works well, but for the HARPS data that we obtained for our sample the Boué et al. (2013) model (as available via the AROME library^{†}^{†}†http://www.astro.up.pt/resources/arome/) should be more appropriate. The model defines line profiles for the CCFs produced by the integrated stellar surface outoftransit, the uncovered stellar surface during transit, and the occulted stellar surface during transit, and assumes them to be even functions. The correction needed to account for the RM effect is calculated through partial differentiation, linearization, and maximization of the likelihood function defined by fitting a Gaussian to the CCF of the uncovered stellar surface. This approach has been tested using simulated data, but has yet to be widely applied to real observations. In this paper, we will therefore compare its results to those from the two other models. To do so, we require values for the width of the Gaussian that is fit to the outoftransit, integrated surface CCF (), and for the width of the spectral lines expected if the star were not rotating (). The latter we set equal to the instrumental profile appropriate to each datum, whilst we use the former as an additional jump parameter for our MCMC algorithm, using the average results given by the HARPS quick reduction pipeline as our initial estimate and applying a prior using that value.
The DT approach was developed by Collier Cameron et al. (2010a) for analysis of hot, rapidly rotating host stars that the RM technique is unable to deal with. It has since been applied to exoplanet hosts with a range of parameters (Collier Cameron et al., 2010b; Brown et al., 2012b; Gandolfi et al., 2012; Bourrier et al., 2015). The alignment of the system is analysed through a comparison of the intransit instrumental line profile with a model of the average outoftransit stellar line profile. This latter model is created by the convolution of a limbdarkened stellar rotation profile, a Gaussian representing the local intrinsic line profile, and a term corresponding to the effect on the line profile of the ‘shadow’ created as the planet transits its host star. This ‘bump’ in the profile is timevariable, and moves through the stellar line profile as the planet moves from transit ingress to transit egress. Its width tells us the width, , of the local line profile, and is a free parameter. Since this width is measured independently, we can disentangle the turbulent velocity distribution of the local profile from the rotational broadening, measuring both and directly. This gives DT an advantage over spectral analysis, as although it is possible to determine the turbulent velocity using the latter method it requires spectra with very high signaltonoise ratio (SNR). For work such as ours it is usually necessary, therefore, to assume a value for .
The path of the bump is dictated by and , and as the planet moves from transit ingress to transit egress its shadow covers regions of the stellar surface with different velocities. This leads to a relation between , , and , which must fit the observed stellar line profile when the local profile and rotational profile are convolved. We thus have two equations for two unknowns ( and , as both and can be determined from the bump’s trajectory), which are therefore well determined. Since and are independently determined using this method, it has the advantage of being able to break degeneracies that can arise between these two parameters in low impact parameter systems (Brown et al., 2012b, e.g.). We note, however, that this breaks down in systems with very slow rotation, i.e. where the uncertainty on is comparable to the rotation velocity.
Another advantage that is often observed with DT is the improved precision on measurements of that it provides, as seen by Bourrier et al. (2015) for the case of the rapidly rotating KOI12 system. This method also has potential as a confirmation method for planetary candidates, as seen with the case of recent case of HATS14 b (Hartman et al., 2015), or conversely as a false positive identifier for difficult to confirm systems.
For all of these models we separate our RV measurements by instrument, and further treat spectroscopic data taken on nights featuring planetary transits as separate data sets. Our Keplerian RV model considers these separated sets of data to be independent. To account for stellar RV noise, an additional m s is added in quadrature to the outoftransit data; this is below the level of precision of the spectrographs used for this work.
2.2 Exploring system architectures
As in our previous work, we explore the possible solutions for each system using a combination of parameter constraints and initial conditions. We have four independent constraints that can be applied.

Apply a Gaussian prior on . This indirectly controls the jump parameters and .

Force the planet’s orbit to be circular, . This indirectly controls the jump parameters and .

Force the barycentric system RV to be constant with time, , neglecting longterm trends that are indicative of third bodies.

Force the stellar radius, , to follow a main sequence relationship with , or use the result from spectral analysis as a prior on .
We consider all 16 possible combinations of these four constraints, analysing each case independently as described above. We discuss these analyses in the following sections. Once all combinations have been examined, we identify the most suitable combination by selecting that which provides the minimal value of the reduced chisquared statistic, . This combination is then reported as the final solution for each system.
3 Results
3.1 Wasp61
WASP61 b orbits a solar metallicity, moderately rotating F7 star, and was initially identified using WASPSouth. Followup observations using TRAPPIST (Jehin et al., 2011), EulerCam (see Lendl et al. 2012 for details of the instrument and data reduction procedure), and CORALIE (Queloz et al., 2000b) showed that the signal was planetary in origin (Hellier et al., 2012). The planet has a circular orbit with a period of d, and has a relatively high density of .
We observed the transit on the night of 2012 December 22 using the HARPS highprecision échelle spectrograph (Mayor et al., 2003) mounted on the 3.6m ESO telescope at La Silla. Fortuitously, we were able to simultaneously observe the same transit photometrically using EulerCam (white light; Fig. 1). We use these new data in conjunction with all of the data presented in the discovery paper (including the original SuperWASP observations) to model the system using our chosen methods.
3.1.1 Hirano model
Trial runs with different combinations of input constraints revealed that the impact parameter of the system is low, . As expected, a degeneracy between and was observed to be present, with the distinct crescent shaped posterior probability distribution covering a wide range of angles and extending out to unphysical values of . Our km s prior restricted the values of the two parameters as expected, but we felt that it was better to allow both to vary normally given our aim of comparing the various RM models. Tests with circular and eccentric orbital solutions showed no evidence for an eccentric orbit, with the Ftest of Lucy & Sweeney (1971) returning a less than percent significance for eccentricity. This is also expected, as our new near and intransit RV measurements do not help to constrain the orbital eccentricity. Relaxing the stellar radius constraint led to insignificant variations in stellar density (which is computed directly from the photometric light curve, and therefore is distinct from the mass and radius calculations), but for some combinations of input constraints the value of the stellar mass varied by . Tests for longterm trends in the barycentric velocity of the system returned results with strongly varying values of both positive and negative , so we set this parameter to zero for our final runs.
The selected solution therefore does not apply a prior on , assumes a circular orbit, neglects the possibility of a longterm trend in RV, and neglects the stellar radius constraint. The fit to our data returns . This particular combination of applied constraints returns an projected spinorbit alignment angle of , an impact parameter of , and a projected rotation velocity of km s, which is in agreement with the spectroscopic value of km s determined from the HARPS spectra using km s, itself derived using the calibration of Doyle et al. (2014). The RM fit produced by the bestfitting parameters is shown in Fig. 2.
3.1.2 Boué model
Similarly to the Hirano model tests, we found no evidence for an eccentric orbit (as expected), no reason to apply a constraint on the stellar radius, and no longterm trend in . The interaction with the prior on was more interesting; with no prior the same degeneracy between and was observed, but while applying the prior restricted the range of rotational velocities explored as expected, it led to a bimodal distribution in . Examination of the posterior probability distribution for the noprior case revealed that this was caused by the Boué model underpredicting compared to the spectroscopic value and the Hirano model, such that the prior from spectral analysis restricted the MCMC algorithm to values within the ‘tails’ of the crescent distribution. Checking posterior distributions for other parameters reveals that the MCMC chain is well converged, and our statistical convergence tests confirm this.
This highlights another degeneracy in the RM modelling problem, in addition to that between and , where orbital configurations with (,) and (,) produce the same ingress and egress velocities and the same chord length (Ohta, Taruya & Suto, 2005; Fabrycky & Winn, 2009). By extension therefore these configurations are indistinguishable when considering the twodimensional problem, and the degeneracy can only be broken by considering the true alignment angle, . This is particularly pernicious in the case of orbits with . Like Fabrycky & Winn, we limit the inclination to the range , which leads to the distribution shown in Fig. 3, with solutions close to .
Ultimately, we adopt the same set of input constraints as for the Hirano model to enable strict comparison between the two models, and acquire final results of km s (consistent with, but lower than the Hirano value as expected from our tests), , and . The resulting RM fit is shown in blue in Fig. 2, which clearly indicates that the two models are fitting the same RM effect in different ways. The Boué fit exhibits steeper ingress and egress gradients, with sharper peaks at larger than the Hirano model, which visually seems to provide a better fit to the RV data although neither model fits the second half of the anomaly particularly well. Comparing their reduced values though reveals that the Boué model gives a slightly poorer fit at , compared to the Hirano model’s .
3.1.3 Doppler tomography
We applied the same set of constraints for our DT analysis as for the other two methods: no prior on ; ; , and no constraint on the stellar radius. The stellar parameters returned were entirely consistent with those from both the Hirano and Boué models. With results of km s, , and , we find no discrepancy between DT analysis and the two other techniques for modelling the RM anomaly. The lefthand panel of Fig. 4 shows the time series of the CCFs, with the prograde signature of the planet barely visible. No sign of stellar activity is visible in the CCF residual map.
Fig. 5 shows the posterior probability distributions for all three analysis methods. Interestingly, in this case tomography seems to provide little improvement in the uncertainties on the alignment angle over the Hirano or Boué models. Instead, the improvement comes in the precision of the measurement, with the uncertainty in the stellar rotation velocity reducing by approximately percent compared to the RM modelling value. This improvement arises due to the different ways in which the different methods treat the spectroscopic data. The Hirano and Boué models are analytic approximations of the behaviour of a Gaussian fit to the composite line profile. This is a valid approach when only the RV data are available, particularly when considering HARPS data as it mimics the calculations performed by the HARPS pipeline. But with the full CCF available, the tomographic method is able to treat the various components of the composite profile explicitly, and can use information from both the timevarying and timeinvariant parts of the CCF directly.
3.2 Wasp62
Like WASP61 b, WASP62 b was discovered through a combination of WASPSouth, EulerCam, and Trappist photometry, in conjunction with spectroscopy from CORALIE (Hellier et al., 2012). The host star is again a solar metallicity, F7type star, and the planet has a circular orbit of period d. WASP62 b is rather inflated ( ) compared to its mass ( ), leading to a much lower density of . Analysis of the HARPS spectra gives km s, with km s from the calibration of Doyle et al. (2014); it is these values that we use for our prior on .
HARPS was used to observe the spectroscopic transit on the night of 2012 October 12. Additional RV measurements were made using the same instruments on 2012 October 15–17 to help constrain the full RV curve. We use the full set of available data to characterize the system, including the weatheraffected EulerCam light curve; as Hellier et al. (2012) note, the MCMC implementation that underlies our analysis accounts for this poorer quality data.
3.2.1 RM modelling
Trial runs to test the effect of applying the four input constraints found that there was no longterm trend in barycentric velocity, and no evidence for an eccentric orbit, with either of the two RM models. Relaxation of the stellar radius constraint led to only minor changes in the reported stellar parameters, with stellar mass and radius being entirely consistent whether the constraint was enforced or not.
Unlike the WASP61 system, the impact parameter was found to be such that no degeneracy was expected between and . This was found to be true for the Hirano model, but our examination of the posterior probability distribution produced using the Boué model showed a long tail in extending out to values that imply very rapid rotation of the host star. In general though, we again find that the Boué model underpredicts compared to the Hirano model  as compared to km s. The Boué model therefore produces larger uncertainties in the value of in order to compensate when trying to fit the RM effect. This effect can be seen in Fig. 8, with the contours for these models being completely distinct. For both models, the alignment angle value remained pleasingly consistent across the different constraint combinations.
The full sets of results, which were produced from runs using no prior on , no constraint on , , and , can be found in Table 5, and show that the larger impact parameter has enabled more stringent limits to be placed on the spinorbit alignment angle than was the case for WASP61. Fig. 6 shows that the two models again produce dissimilarly shaped bestfitting RM models; as with WASP61, the Boué model has steeper ingress and egress velocity gradients, and sharper peaks. However, the angles produced by the two models are entirely consistent, with the Hirano model finding and the Boué model .
3.2.2 Doppler tomography
Using the same set of input constraints as for our RM modelling, we again carried out DT analysis of the system, finding an alignment angle of . The planetary signature of WASP62 is much stronger than that of WASP61, and can be seen far more clearly in the CCF time series map (Fig. 7). For this system, the tomographic method has improved the uncertainties on by roughly a factor of compared to the Boué model, but again provides little improvement over the precision afforded by the Hirano model. All three results are consistent with alignment according to the criterion of Triaud et al. (2010), a conclusion which is supported by the trajectory of the planetary signal in Fig. 7. However, this is based on an ad hoc criterion that reflects the typical uncertainty on measurements at the time that it was formulated. Work since 2010 has improved the typical uncertainty, such that this criterion is no longer really applicable. We therefore classify WASP62 as slightly misaligned; the resolution of the planet trajectory in our Doppler map is insufficient to distinguish this from a truly aligned orbit.
As with WASP61, it is the treatment of by the three models that is interesting here. We have already noted that the Boué model returns lower values than the Hirano model, but the DT result of falls between the two whilst being consistent with neither thanks to the small error bars on all three estimates (see Fig. 8). None of the values that we find are consistent with the spectroscopic value of km s derived from the HARPS spectra. We also note that the uncertainties in using this method are smaller than for either the Hirano or Boué models (see Table 5).
3.3 Wasp71
Smith et al. (2013) presented the discovery of WASP71 b using photometry from WASPN, WASPSouth and TRAPPIST, along with spectroscopy from CORALIE that included observations during transit made simultaneously with the TRAPPIST observations. The host star was found to be an evolved F8type, and significantly larger and more massive than the Sun, whilst the planet was found to be inflated compared to the predictions of Bodenheimer, Laughlin & Lin (2003, e.g.), and to have a circular orbit with a period of d. The spectroscopic transit observations made using CORALIE enabled Smith et al. to measure the projected spinorbit alignment angle of the system. They found that the system was aligned, with , and rapidly rotating at km s (calculated assuming km s following Doyle et al. 2014).
We obtained additional spectroscopic data on the night of 2012 October 26, observing a complete transit, with further observations made on 2012 October 23 and 25. We combine these with the discovery photometry and spectroscopy to model the system. We do not, however, include the spectroscopic transit used by Smith et al. to measure , for two reasons. The first is that we wish to obtain an independent measurement of the spinorbit alignment. The second is that, as noted by Boué et al. (2013), different instruments can produce different signals from the same measurement owing to their different analysis routines, and therefore RM data sets from different instruments should not be combined. Analysis of our new HARPS spectra gives km s and km s.
3.3.1 RM modelling
Using the Ftest of Lucy & Sweeney (1971) we found no indication of significant eccentricity in the system, in agreement with Smith et al. (2013), and therefore set in our final analysis. We also found no consistent evidence that there is a longterm trend in barycentric velocity, so set .
The interaction between , , and the stellar radius constraint is an interesting one for this system. Relaxing the constraint on causes the stellar radius to decrease by percent, with the stellar mass increasing by approximately percent. Relaxing the constraint also leads to a significant, approximately tenfold rise in the impact parameter from to , with corresponding effect on the result for , which with the radius constraint active is almost unphysically large owing to the degeneracy that arises with both and . Smith et al. (2013) report an impact parameter of , so we chose not impose the stellar radius constraint to allow the impact parameter to fit to what appears to be the more natural value. This does mean that we find a larger, less dense planet than the Smith et al. (2013) result. We also note that applying the stellar radius constraint returns a stellar effective temperature which is K hotter than previous spectroscopic values, whilst neglecting the constraint gives a temperature more consistent with previous analyses.
Once again, we find that the two different methods return similar results for the alignment angle and stellar parameters, but that the Boué model gives a more slowly rotating star than is suggested by the Hirano model (see Table 5). Fig. 9 shows the bestfitting models produced by both methods, with the Hirano model (the dashed, red line) having a shallower peak during the first half of the anomaly. There is substantial scatter in the RV measurements during this period however, and the Hirano model appears to better fit the second half of the anomaly, where there is less scatter in the radial velocities.
The final solutions that we report were taken from runs with no prior on , , , and no constraint applied to the stellar radius.
3.3.2 Doppler tomography
Whereas for the two previous systems the tomographic analysis supported the Hirano model with regards to the projected rotation velocity of the host star, for WASP71 it is the Boué model with which DT agrees (see Fig. 10), although the value of km s that we find is significantly lower than the spectroscopic value. The other parameter values that are found through DT are substantially different to either set of RM results. The impact parameter is lower, leading to a lower value of that is consistent with (see Table 5). Particularly interesting though is the difference in the physical stellar parameters found by this method, which imply a smaller star. As implied by our result of , Fig. 11 shows that the system is well characterizedaligned, in agreement with the result from Smith et al. (2013), although we are not able to significantly improve on the precision that they report.
3.4 Wasp76
WASP76 A (West et al., 2016) is another F7type planethosting star, but is rotating significantly more slowly than either WASP61 or WASP62. The planet is substantially bloated, with a density of only , and orbits its host every days in a circular orbit. It was discovered and characterized using data from WASPSouth, TRAPPIST, EulerCam, SOPHIE (Bouchy et al., 2009; Perruchot et al., 2011), and CORALIE.
HARPS was used to observe the transit taking place on 2012 November 11, and to make additional measurements on 2012 November 12–14. We combined these measurements with the discovery paper’s photometry for our analysis, excluding two spectra that were obtained at twilight. Spectral analysis of the new spectra returned km s using the calibration of Doyle et al. (2014), leading to km s. We use this for our prior on rotation velocity. The SNR of the spectra are relatively poor however, so we increased the lengths of our MCMC phases to (minimum ) for burnin, and for the production phase, leading to a final chain length of for the concatenated chain. We also approached the analysis of this system in a different manner to the other systems in our sample.
3.4.1 Hirano model
We began by testing the effects of constraints , , and (see Section 2.2). We found no evidence for a longterm trend in barycentric RV, so adopted the constraint for our final solution. We also set after finding no evidence for a significantly eccentric orbit. Despite the large number of photometric light curves available for the system, we found that applying the stellar radius constraint led to increases in both and . We attribute this to a combination of poor photometric coverage of the transit ingress, and significant scatter in some of the light curves (see fig. 1 of West et al. 2016). Despite the varying stellar parameters there was no compelling reason to apply the constraint, and we therefore chose not to do so for the next phase of our analysis.
Exploring
Initial exploratory runs with , , no radius constraint, and no prior on consistently found a small impact parameter of , in agreement with the discovery paper value of but poorly constrained. These runs produced the expected degeneracy between and , and the associated crescentshaped posterior probability distributions (see Fig. 14). The distribution shows a slight preference for positive , but the uncertainty on the value was large. Furthermore, our Geweke and Gelman–Rubin tests implied that the MCMC chains were poorly converged.
We thus applied constraint , a prior on . The addition of this constraint leads to a bimodal distribution in , with both minima being tightly constrained (see Fig. 14). Further investigation revealed that individual chains were split roughly between the positive and negative minima in space, dependent on the chain’s exploration of parameter space during the burnin phase; this naturally led to poor convergence when analysing the concatenated MCMC chain, but inspection of trace plots, autocorrelation data, running means, and statistics from the Geweke test showed that each individual chain was well converged. We thus ran additional chains to collect five that favoured the positive minimum and five that favoured the negative minimum, and tested the convergence of the two minima. Both were found to be well converged. We also carried out tests whereby the chain was started at ; the results matched our expectations, with each chain remaining in the associated positive / negative minimum and being well converged.
As with some of our modelling of the WASP61 (see Section 3.1.2), this is an example of the degeneracy inherent in the RM problem, whereby solutions with (, ) are indistinguishable in terms of fitting the data from solutions with (, ).
Interestingly, chains that explored the positive minimum consistently returned a higher value of the impact parameter than those chains that explored the negative minimum. However, in the case of the positive minimum the median impact parameter of was consistent only with the lower end of the impact parameter given in the discovery paper, , and in the case of the negative minimum the value of did not agree with discovery paper at all.
Impact parameter
We therefore elected to explore the option of applying an additional constraint on the system, this time on the impact parameter using the value from the discovery paper as a prior. We tested the application of this prior to cases both with and without a prior on .
When we applied the prior on but not the prior on , we found very similar results to those obtained in the corresponding case without the impact parameter prior, albeit with much improved convergence of our MCMC chains. The concatenated chain showed a preference for the minimum, with sizeable uncertainty on the result; the major difference with that earlier example was the more tightly constrained impact parameter distribution. Testing chains with initial alignment of showed completely consistent results with the free case, with both cases favouring the positive minimum, albeit with substantial uncertainty on .
When priors on both the impact parameter and were both applied to the free case, the MCMC chains were forced into the positive minimum. The results from the concatenated chain give an impact parameter in agreement with the discovery paper’s value, and in addition provide a more precise determination of than any of the other combinations of constraints. The gave solutions in the corresponding minima, but it is notable that the impact parameter for the negative case is significantly lower than the value expected from the discovery paper. This suggests that the positive solution should be favoured.
Results from these analyses are shown in Table 3.
3.4.2 Boué model
We investigate the system using the Boué model, following the same methodology outline for the Hirano model. We again adopt a constraint of zero drift in the barycentric velocity owing to lack of evidence to the contrary. Applying the stellar radius constraint led to increases in , , , , , and , but in several cases these parameters were unphysical. There was also no substantial improvement in fit when applying the radius constraint, and we therefore elected not to do so. We also adopted a circular solution as there was no evidence for a significantly eccentric orbit. In this we match our choice of constraints for the Hirano model, which is encouraging as it again shows that the two models are broadly consistent in their exploration of parameter space.
Initial tests without a prior on also produced results consistent with those found using the Hirano model, including the crescentshaped degeneracy between and with a preference for positive , though the stellar rotation velocity was found to be even slower at km s. When we applied the prior, we found that although convergence statistics were greatly improved, and the rotation velocity now agreed with our expectations from spectral analysis, the impact parameter was significantly lower than expected, and a bimodal distribution in was obtained. When we forced , the convergence statistics were again improved and the chains explored the expected minimum, but like the Hirano model tests the impact parameter remained lower than anticipated.
Applying a prior on the impact parameter, in the absence of the prior, showed that the chains favoured the positive minimum, but with sizeable uncertainty on and a slow rotation velocity. When we apply priors on both impact parameter and , we found results consistent with the Hirano model equivalents.
Model  

prior?  prior?  /  /km s  /  /  
Hirano  No  No  
Yes  No  
Yes  No  
Yes  No  
No  Yes  
No  Yes  
No  Yes  
Yes  Yes  
Yes  Yes  
Yes  Yes  
Boué  No  No  
Yes  No  
Yes  No  
Yes  No  
No  Yes  
No  Yes  
No  Yes  
Yes  Yes  
Yes  Yes  
Yes  Yes  
Tomography  No  No  
Yes  No  
Yes  No  
Yes  No  
No  Yes  
No  Yes  
No  Yes  
Yes  Yes  
Yes  Yes  
Yes  Yes 
3.4.3 Doppler tomography
We adopted constraints of and , but left the stellar mass and radius freely varying, in order to be consistent with our analyses using the Hirano and Boué models.
As anticipated, using the DT method substantially reduced the degeneracy between and , though it did not, for this system, remove it completely (see Fig. 14). DT again favoured the positive minimum, though more strongly than the other two models, and again returned a more slowly rotating star than anticipated. Adding a prior on , however, produced different behaviour than shown previously. With DT, adding a prior on forced the chains into the negative minimum, with no bimodal distribution observed, though once again the impact parameter strongly disagreed with the value from the discovery paper and implied a central transit.
If we impose a prior on the impact parameter using the value from the discovery paper, then in the absence of a prior on we again find that the chains favour the positive minimum in , irrespective of the value of . We also find a faster value of than was returned by either the Hirano or Boué models for the same combination of priors, though the values are consistent to . If we add the prior on then the results remain consistent, but the uncertainties are reduced in magnitude, particularly for which also moves closer towards a value that implies a polar orbit.
3.4.4 A possible polar orbit?
Consecutive analyses of WASP76 have gradually reduced our assessment of the stellar rotation velocity. Spectral analysis of the CORALIE data by West et al. (2016) gave km s, while analysis of our new HARPS spectra returned a value of km s. In the absence of a prior on rotation, our modelling of the RM effect using any of the three methods gives a value significantly slower than this, at km s. Yet inspection of both the CORALIE and HARPS spectra reveals visible rotation (see Fig. 15), and the full width at halfmaximum (FWHM) of the spectra are greater than for stars with similar () colour, such as WASP20. This would suggest that the star is indeed oriented close to edgeon, rather than the poleon orientation suggested by the slow rotation velocity returned by our MCMC chains.
We therefore consider the possibility that the orbit is oriented at close to , with a transit chord such that the path of the planet is almost parallel to the stellar rotation axis. This solution is consistent with the path of the planetary ‘bump’ through the stellar line profile in Fig. 13, which shows little movement in velocity space. To explore this possible system configuration, we carried out additional analyses both with and without a prior on , this time forcing the MCMC chain to adopt throughout. Note that these analyses used the constraints of and , and applied no constraints on the stellar mass or radius, as before. We present these results in Table 4 and Fig. 12.
Model  

prior?  prior?  /  /km s  /  /  
Hirano  Yes  No  (fixed)  
(fixed)  
Yes  (fixed)  
(fixed)  
No  No  (fixed)  
(fixed)  
Yes  (fixed)  
(fixed)  
Boué  Yes  No  (fixed)  
(fixed)  
Yes  (fixed)  
(fixed)  
No  No  (fixed)  
(fixed)  
Yes  (fixed)  
(fixed)  
Tomography  Yes  No  (fixed)  
(fixed)  
Yes  (fixed)  
(fixed)  
No  No  (fixed)  
(fixed)  
Yes  (fixed)  
(fixed) 
We found that when applying the Boué model, the MCMC chains took approximately twice as long to converge as when applying the Hirano model. The source of this difficulty with convergence is uncertain, but seems to be related to the ratio between and . Large steps in the latter that explore rapidly rotating solutions lead to unphysical values of this ratio, such that the step fails. This restricts the set of possible solutions to a more limited area of parameter space, such that a larger percentage of possible steps lead to poor solutions, and thus convergence of the chain proceeds more slowly.
The three analysis techniques generally produced consistent results. In the majority of cases, we found that the value for returned by the chains was in agreement with the discovery paper; the exceptions to this were the cases with and the prior only. With the rotation prior inactive, the only case to be consistently in agreement with predictions across all three techniques was the case with prior also inactive, and ; in the other cases, the stellar rotation was generally slower than the spectral analysis result (as noted in previous sections).
These results lend some small support to the hypothesis of a polar orbit, as we note that in the case of neither prior being applied the results were consistent with both spectral analysis and the discovery paper.
3.4.5 A poorly constrained system?
WASP76 seems to represent a similar case to WASP1 (Albrecht et al., 2011): a low impact parameter, combined with a poor SNR for the RM effect, leading to a weak detection. Here though, we find a threeway degeneracy between , , and that can only be broken through the application of appropriate Gaussian priors.
Our results show tentative support for a strongly misaligned orbit. Applying a prior on stellar rotation or on the impact parameter produces results that suggest strong misalignment, particularly when using the DT method. Forcing the system to adopt a polar orbit, or using a polar orbit as the initial condition, reveals that this is a plausible option for the system’s configuration, though there is still substantial ambiguity in the precise orientation of the planet’s orbit.
Although we have reported and discussed results for cases both with and without the various combinations of these two priors, in Section 4 we focus on the case with a prior on the impact parameter, but without a prior on . This combination maximizes the relevance of our crosssystem comparison by allowing the different models to evaluate stellar rotation freely, while ensuring that the other system parameters are truly representative and derived from fully converged chains; for the WASP76 system this necessitates the prior on . We caution readers, however, that we cannot constrain the obliquity of the planet’s orbit beyond the general statement that it is likely strongly misaligned in the positive , prograde direction.
3.5 Wasp78
WASP78 b (Smalley et al., 2012) is largeradius, lowdensity hot Jupiter that orbits its host star every d in a circular orbit. From the CORALIE RV data in the discovery paper, the orbit appears to be circular, and such a solution is adopted therein, but the authors state that additional RV data is required to pindown the eccentricity. Spectral analysis shows that the host star is of spectral type F8, with K, km s, and km s.
We obtained new photometry using EulerCam, in both the and bands, of transits on 2012 November 2 and 26, respectively, which we combine with the discovery photometry from WASPSouth and TRAPPIST. The band light curve shows clear modulation, almost certainly as a result of the presence of starspots or related stellar activity features, as this modulation is not replicated in the band light curve (Fig. 17). Our HARPS spectroscopic transit observations were made simultaneously with the band EulerCam observations, and used in conjunction with spectroscopy from CORALIE. Analysis of the HARPS spectra provides km s, and km s.
3.5.1 RM modelling
As with the other systems in our sample, we find no evidence of a longterm barycentric velocity trend, and set . We also find little difference in behaviour or results between the two models; the following discussion is applicable to both.
When allowing eccentricity to float we find that the MCMC algorithm returns two general sets of solutions, with the Hirano and Boué models giving comparable estimates as expected. The imposition of the stellar radius constraint leads to , but forces the stellar radius to be approximately half the value presented in Smalley et al. (2012). Furthermore, examination of Fig. 18 shows that high eccentricity orbits are a very poor fit to the RV data. We do not consider these solutions plausible. When the stellar radius constraint is removed, the free eccentricity fit returns . Testing these small eccentricity solutions using the method of Lucy & Sweeney (1971) reveals that they are not significant. Comparing eccentric and circular solutions with no other constraints applied, we find that variations in physical parameters are within the uncertainties. We thus choose to force a circular orbit as it provides the most plausible solution; the phase coverage of the HARPS data is insufficient to truly constrain any small eccentricity that might be present, and the preexisting CORALIE RV data cannot provide a firm conclusion regarding circularity (or otherwise). We also note that hot Jupiters with M are generally observed to have zero (or small) eccentricity (Anderson et al., 2012).
Application of the stellar radius constraint to a circular orbit increases the reported stellar mass, and also leads to an increased value of and a lower impact parameter. Our results in this scenario agree very well with those from Smalley et al. (2012), as expected; the host star is both too massive and too large, and inconsistent with both spectral analysis of both the existing CORALIE spectra and our new HARPS data. Conversely, if the planet’s orbit is allowed to be eccentric (against the evidence), then adding the stellar radius constraint decreases the mass of both bodies, and significantly decreases their radii.
We find a moderate impact parameter of , with no correlation or degeneracy between and , so do not apply a prior using the spectroscopic . Our chosen solution is therefore a circular orbit with no longterm velocity trend, and no application of the constraints on or . The full set of results is shown in Table 5, but the alignment angles found by the Hirano and Boué models are consistent with each other, and with , implying a wellaligned orbit. We do however find that the Boué model produces a slower rotation velocity, as we found for all of the previous systems, in this case inconsistent with the spectroscopic value. This cannot, however, be as a result of a poorly constrained impact parameters; the results in Table 5 show that the impact parameter returned by the Boué model is in agreement with the results from other models, and has similar uncertainties.
The form of the bestfittingting models are very similar in Fig. 19, as they were for WASP76 (the other system with a low signaltonoise anomaly), although in this case there are far fewer data during the transit.
Our selected solutions are taken from runs for which no prior was applied to , , , and no constraint was applied to the stellar radius.
3.5.2 Doppler tomography
We apply the same set of constraints to our tomographic MCMC analysis (Fig. 20) as for the RM modelling runs. We find that all parameters are in agreement with the results from the Hirano model and the case in which no RM modelling is carried out, and that all parameters except are in agreement with the Boué model results, though we note that the disagreement there originates with the Boué model rather than DT. Unlike for some of the other systems studied herein, we find a roughly factor of improvement in the precision of our alignment angle measurement over the RM modelling methods when using tomographic analysis (see Fig. 21 and Table 5).
Acquiring simultaneous photometry and spectroscopy provides us with a means to crosscheck for stellar activity signatures, as any stellar activity which affects the light curve should be visible in the CCF time series map, Fig. 20. Fig. 17 shows evidence of the presence of star spots in the band light curve, but the band light curve that was acquired simultaneously with our HARPS spectroscopy shows only hints of similar activity. Examination of the righthand panel of Fig. 20 is inconclusive.
3.6 Wasp79
Smalley et al. (2012) found WASP79 b to have both a low density and a large radius, and to orbit its F5type host star in a circular orbit with a period of d. But they were unable to fully characterize the shape and duration of the transit signature owing to a lack of highprecision followup photometry, having access to only a single, partial transit light curve from TRAPPIST. This limited the accuracy of the physical and orbital parameters that they were able to obtain. Addison et al. (2013) measured the spinorbit alignment angle of the system using UCLES, modelling the RM effect during transit to obtain , indicating significant misalignment.
We obtained band EulerCam observations of the transits on the nights of 2012 November 11 and December 4 (see Fig. 22), as well as spectroscopic observations of the transit on 2012 November 13 using HARPS. We also used the data from WASPSouth, TRAPPIST, and CORALIE that was presented in Smalley et al.. Spectral analysis of the HARPS spectra provides values of km s and km s that we use as our prior on the rotation velocity. We note though that the value of is extrapolated owing to the star’s effective temperature being outwith the range used for the calibration of Doyle et al. (2014). If is fixed to , then we obtain km , which effectively places an upper limit on the stellar rotation velocity of km s.
3.6.1 RM modelling
We find no evidence for an eccentric orbit or for a longterm barycentric velocity trend, and find no difference between the stellar parameters obtained by MCMC when run with and without the stellar radius constraint. We also find no reason to apply a prior on ; the impact parameter is sufficiently large for there to be no discrepancy between and , although we do once more see a more slowly rotating star with the Boué model than with the Hirano model. In this case the Boué model is consistent with the spectroscopic rotation velocity, whilst the Hirano model seems to be overpredicting the speed of the stellar surface, even when we consider the upper limit set by the zero macroturbulence case. Of interest are the different forms of the [] posterior probability distributions for the two models; that for the Hirano model shows a correlation between the two parameters and a more triangular shape, whilst the distribution for the Boué model is more elliptical in shape.
Our selected solutions are taken from runs for which no prior was applied to , , , and no constraint was applied to the stellar radius. Both models return strongly misaligned results for the system at very high precision: for the Hirano model, and for the Boué model, in agreement with the result of Addison et al. (2013). The reason for such strong constraints are readily apparent from Fig. 23. The RM anomaly is strongly asymmetrical and comprises only positive deviation from the outoftransit velocity. This indicates that only one half of the star, the approaching hemisphere, is being traversed by the planet, and in combination with the angle that we find suggests a nearpolar orbit for the planet.
3.6.2 Doppler tomography
Our tomographic analysis is similarly well constrained, and returns an alignment angle of , consistent with the results from both the Boué and Hirano models. What is notable however is that the uncertainties have been reduced by a factor of over the already impressive results that we obtained with those models, and are an order of magnitude better than those obtained by Addison et al. (2013). for our DT analysis falls between those of the two RM models, and is faster than the spectroscopic result, though it is in agreement with the upper limit set by the zero macroturbulence case (see Fig. 24 and Table 5).
Fig. 25 again shows the strongly asymmetric signature of a misaligned polar orbit. The planetary trajectory is unlike that for any of the other five systems, being confined to one side of the stellar spectral line and moving from right to left (although any movement through the line is slight at best). This implies either a polar orbit, or a slightly retrograde one, as suggested by the value of that we obtain.
4 Discussion
We have obtained measurements of , , and using three different analyses of the RM effect. We show these results in Table 5, along with other relevant parameters as determined from either spectral analysis of our new HARPS data (), modelling of the available photometric and RV data (, ), or isochronal fitting (, see Section 4.2). In the following discussion we will consider our results in the context of existing literature data on spinorbit misalignment, treating the values obtained through the DT technique as definitive. We provide a more comprehensive list of parameters in Table 7 in the appendix.
System  Model  /km s  /  /  /  /  /Gyr  

WASP61  Hirano  
Boué  
Tomography    
No RM fitting      
WASP62  Hirano  
Boué  
Tomography    
No RM fitting      
WASP71  Hirano  
Boué  
Tomography    
No RM fitting      
WASP76  Hirano  
Boué  
Tomography    
No RM fitting      
WASP78  Hirano  
Boué  
Tomography    
No RM fitting      
WASP79  Hirano  
Boué 