{'fraction_non_alphanumeric': 0.027024021352313167, 'fraction_numerical': 0.019943653618030844, 'mean_word_length': 4.50214154599225, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Line-profile tomography of exoplanet transits I: The Doppler shadow of HD 189733b

We present a direct method for isolating the component of the starlight blocked by a planet as it transits its host star, and apply it to spectra of the bright transiting planet HD 189733b. We model the global shape of the stellar cross-correlation function as the convolution of a limb-darkened rotation profile and a gaussian representing the Doppler core of the average photospheric line profile. The light blocked by the planet during the transit is a gaussian of the same intrinsic width, whose trajectory across the line profile yields a precise measure of the misalignment angle and an independent measure of v sin I. We show that even when v sin I is less than the width of the intrinsic line profile, the travelling Doppler"shadow"cast by the planet creates an identifiable distortion in the line profiles which is amenable to direct modelling. Direct measurement of the trajectory of the missing starlight yields self-consistent measures of the projected stellar rotation rate, the intrinsic width of the mean local photospheric line profile, the projected spin-orbit misalignment angle, and the system's centre-of-mass velocity. Combined with the photometric rotation period, the results give a geometrical measure of the stellar radius which agrees closely with values obtained from high-precision transit photometry if a small amount of differential rotation is present in the stellar photosphere.


INTRODUCTION
The recent discoveries of strong projected spin-orbit misalignments in the transiting exoplanets XO-3b (Hébrard et al. 2008, Winn et al. 2009a, HD 80606b (Winn et al. 2009b, Pont et al. 2009) and WASP-14b ), and retrograde orbital motion in WASP-17b (Anderson et al. across the stellar disc, the projected misalignment angle λ between the orbital axis and the stellar spin axis, and the relative sizes of the star and planet. Detailed semi-analytic formulations for determining the anomalous departure of the line centroid from the stellar reflex orbit have been developed by Ohta, Taruya, & Suto (2005), Giménez (2006) and Hirano et al. (2009). The process of measuring and interpreting this reflex velocity shift is not, however, straightforward. If the stellar rotation profile is even partially resolved by the spectrograph, the missing light within the planet's silhouette introduces an asymmetry into the line spread function. For instruments such as SOPHIE and HARPS, radial velocities are derived using a gaussian fit to a cross-correlation function (CCF) computed from the stellar spectrum and a weighted line mask (Baranne et al. 1996;Pepe et al. 2002). Such methods work extremely well for the symmetric and time-invariant profile shapes encountered outside transit. When applied to profiles with a time-varying degree of intrinsic asymmetry, however, they yield velocity measures that depart systematically from the velocity centroid of the visible starlight. Winn et al. (2005) overcame this difficulty for their iodine-cell observations of HD 209458b using a direct modelling approach which they have employed in subsequent papers. They built models of both the unobscured stellar profile and the light blocked by the planet into their analysis of the line-spread function, deriving semi-empirical corrections to the model velocities to enable meaningful comparison with the data. Nonetheless, their study of the RM effect in HD 189733 (Winn et al. 2006) shows a clear pattern of correlated radial-velocity residuals during the transit. In a more recent HARPS study of the RM effect during a transit of HD 189733b, Triaud et al. (2009) noted an almost identical systematic pattern of correlated residuals between the radial velocities of the gaussian fits to the CCFs, and the line centroid velocities of the best-fitting model computed using the formulation of Giménez (2006).
In developing an analytic model of this velocity anomaly for cross-correlation spectra, Hirano et al. (2009) have shown that such errors are exacerbated as v sin I increases and the asymmetry becomes more extreme. If uncorrected, the anomaly can lead to significant over-estimation of the value of v sin I, particularly for rapidly-rotating planet-host stars. The motivation for the present study is to develop a comprehensive model of the changes in CCF morphology that occur before, during and after a transit, in terms of the projected stellar equatorial rotation speed, the spin-orbit misalignment angle, the intrinsic width of the combined stellar CCF and instrumental broadening function, and the stellar limb-darkening coefficient. In principle this approach is similar to that of Winn et al. (2005), but we use it to decompose the observed CCF into its various components and track them directly. A similar methodology has already been used to model the RM effect in the binary stars V1143 Cyg (Albrecht et al. 2007) andDI Her Albrecht et al. (2009). To illustrate the effectiveness of this approach as an alternative to gaussian-fitting in the study of transiting exoplanets, we re-analyse the same HARPS observations of HD 189733b described by Triaud et al. (2009). We show that the spectral signature of the light blocked by the planet is clearly discernible in the residuals when a model of the unobscured starlight is subtracted from the data, and track it directly to derive new measurements of the projected stellar rotation velocity and spin-orbit misalignment angle.

OBSERVATIONS AND DATA ANALYSIS
Densely-sampled spectral time-series of HD 189733 were obtained with the HARPS instrument on the ESO 3.6-m telescope at La Silla on the nights of 2006 September 8 and 2007 August 29. The observations are described in detail by Triaud et al. (2009). The 2006 September data were taken with a cadence of 10.5 minutes, while the 2007 data were sampled every 5.5 minutes. The signal-to-noise ratio (SNR) in the continuum of each spectrum was typically 230 in echelle orders 61 through 70 of all spectra in the 2006 sequence, and 130 for the 2007 data which employed shorter exposures.
The spectra were extracted and numerical crosscorrelation functions (CCF) were computed using the most recent version of the HARPS data-reduction software (Mayor et al. 2009). This involves order-by-order crosscorrelation of the extracted stellar spectrum with a weighted mask function, as described by Baranne et al. (1996) and Pepe et al. (2002). The mask used for these observations was derived from the line list for a star of spectral type K5. The SNR in the pseudo-continuum to either side of the cross-correlation peak was measured to be ∼ 6300 for a typical spectrum in the 2006 series, and ∼ 3900 in the 2007 sequence. The cross-correlation process thus yields a multiplex gain factor of approximately 29 in SNR, for both data sets, allowing the SNR in any individual CCF to be estimated directly from the SNR of the original spectrum.

Stellar CCF profile model
The CCF of a slowly-rotating star exhibits a flat continuum and a single dip, which can be fitted to high precision with a gaussian profile. We assume that the local line profile at any point on the surface of a limb-darkened, rotating star takes the form of a gaussian, where the dimensionless velocity x = v/v sin I and gaussian sigma s = σ/v sin I are expressed in units of the projected stellar equatorial rotation speed v sin I. We assume a standard linear limb darkening model with limb-darkening coefficient u: using the standard convention that the direction cosine of the local surface normal relative to the line of sight is µ = cos θ, with θ = 0 at the centre of the stellar disc. Assuming solid-body rotation, the limb-darkened rotation profile is obtained by substituting µ = p 1 − x 2 − y 2 and integrating Eq. 2 between limits y = ± √ 1 − x 2 , where x = v/v sin I as before. The resulting expression is normalised by integrating over the range −1 < x < 1 to obtain (Gray 1976) These expressions are normalised such that The observed rotation profile is the convolution of the two functions, which is computed via a straightforward numerical integration. Since the CCFs are computed in the velocity frame of the solar-system barycentre, the shift of the model CCF to match a Keplerian orbit solution is achieved by defining xij for a pixel i with velocity vij in the barycentric frame to be where νj is the true anomaly at the time of the jth observation and all other orbital elements K, e,ω, γ have their usual meanings. The model of the jth observation of the stellar CCF is thus h(xij).

Missing-light model during transit
At any moment, the position of the planet on the plane of the sky is given by where r is the instantaneous distance of the planet from the star, and i is the inclination of the orbital angular momentum vector to the line of sight. The z axis is parallel to the projected direction of the orbital axis ( Fig. 1). If xp, zp and r are expressed in units of the stellar radius, the fraction of The model stellar profile h ij has been subtracted, leaving the bright travelling signature of the starlight blocked by the planet during the transit. This is seen superimposed on a fixed-pattern residual spectrum α i which remains stable during each night but changes subtly from one night to the next.The horizontal dotted line gives the phase of mid transit. Crosses denote times of first, second, third and fourth contact. The stellar radial velocity is plotted as a dotted curve flanked by dashed lines at ±v sin I.
the starlight blocked by the planet during the total part of a transit is where µ = p 1 − x 2 p − z 2 p . The radial velocity of the missing starlight relative to the star's instantaneous velocity depends on the projected distance of the planet from the stellar rotation axis. Following Winn et al. (2005) we define λ = φspin − φ orbit where φ denotes position angle in the plane of the sky. The scaled best-fit model h ij of the stellar profile, the travelling planet signature β g(x ij − up) and the fixed sidelobe pattern α i have all been subtracted from the data as in Eq. 16, leaving only noise.
Rotating the coordinate frame in the plane of the sky, we obtain the perpendicular distance of the planet from the stellar rotation axis in units of the stellar radius: which in turn yields the barycentric velocity of the missing light, vp = upv sin I. The combined model of the stellar CCF and the missing starlight during the total part of the transit is thus Analytic expressions for computing β and up during the partial phases of the transit are given by Ohta, Taruya, & Suto (2005) and Giménez (2006). The term β g(xij −up) represents the travelling gaussian "bump" produced by the blocked starlight. The appearance of this residual bump after subtraction of the model stellar profile from the observed CCFs is illustrated in greyscale form in Fig. 2. The coefficient β depends on the ratio Rp/R * and the limb-darkening coefficient u. The width of the gaussian g gives an additional constraint on the non-rotating local CCF width s, given the reasonable assumption that Rp << R * . The trajectory of the planet signature through the line profile is described by up, which in turn depends on v sin I and λ. These separate constraints on s, u and v sin I are valuable, since they help to break any degeneracy in the description of the shape of the stellar CCF, particularly when v sin I is small compared to the intrinsic linewidth.

Model fitting and error scaling
The data and model are both orthogonalised by subtracting their optimal mean values: using inverse-variance weights wij = 1/σ 2 ij . The multiplicative constantÂ by which the model must be scaled to give an optimal fit to the data is then given by the expression The CCF obtained by cross-correlating a late-type stellar spectrum with a line mask inevitably contains weak sidelobes, caused by random alignments between stellar lines and mask lines at arbitrary shifts. The cross-correlation dip therefore appears superimposed on a background pseudocontinuum which exhibits weak fixed-pattern residuals. These appear as an alternating pattern of light and dark vertical stripes in Fig. 2. Imperfect modelling and removal of the spectrograph blaze function prior to cross-correlation can also produce a weak asymmetry in the CCF profile, which would leave a similar pattern in the residuals. At this stage in the calculation, the sidelobe pattern is removed by computing and subtracting the optimal average of the residual spectra: The sidelobe pattern can also be influenced by moonlight and telluric absorption levels, and potentially by weak, timevariable asymmetries in the CCF caused by starspots. This produces subtle changes in the sidelobe pattern from one night to the next. Little change is seen within a single night. The method is therefore most effective when the spectrumaveraging summations are performed separately for different nights. The residuals following the subtraction of the sidelobes and the signature of the light blocked by the planet are shown in Fig. 3.
The goodness of the model fit to the observed CCF is then measured via the statistic Care is needed when estimating the error bars on both the data and the fitted parameters. The spectral resolving power of HARPS is R = λ/∆λ = 120000, yielding a velocity resolution of 2.5 km s −1 . The CCFs produced by the HARPS data reduction pipeline are oversampled, being binned in velocity increments of 0.25 km s −1 . The errors on adjacent points in the CCF are therefore strongly correlated. In order to ensure that all data points contributing to χ 2 are statistically independent, we bin the data and the fitted model in velocity, by a factor ten before computing χ 2 . Initially we estimate the relative uncertainties σij of the CCF data in each spectrum from the SNR of the original spectrum using the multiplex gain factor described in Section 2 above.

Markov-chain Monte Carlo parameter fitting
Given a sequence of CCFs derived from a set of spectra densely sampled during transit and distributed around the rest of the orbit, the overall χ 2 statistic is therefore a function of the orbital parameters K, e, ω and γ; the orbital period P and the epoch Tp of periastron; the parameters v sin I, s and u describing the shape of the rotationallybroadened stellar CCF and the continuum tilt; and the error scaling factor f . The limited phase range of the HARPS observations does not permit as reliable a determination of the full set of orbital parameters as the combined SOPHIE and HARPS data used by Triaud et al. (2009). We therefore chose not to fit the full radial-velocity orbit and the photometric transit parameters simultaneously with the spectroscopic transit parameters. Instead, we fixed the values of all other model parameters at the values given by Triaud et al. (2009) in their recent combined photometric and spectroscopic analysis of the same system. These values are summarised and compared with their values from earlier studies in Table 1. Only the system centre-of-mass velocity, γ, was re-determined directly from the data. This was done because the star is known to exhibit starspot activity, so small nightto-night changes in γ are to be expected as the star rotates. Although spot activity could in principle affect the apparent value of K too, the latter is well-determined from RV observations collected over many orbits. Any degradation of the fit by starpot-induced changes in the apparent radial acceleration during the transit is expected to be small. We determined the model parameters v sin I, λ, s and γ and their joint posterior probability distribution using a Markov-chain Monte Carlo (MCMC) approach.
Strongly-correlated pairs of parameters are best avoided in MCMC, since they lead to slow convergence and undesirably strong correlations between successive steps in the chain (Tegmark et al. 2004). We expect the parameters s and v sin I, which together determine the width of the stellar CCF outside transit, to be correlated to some extent, since their quadrature sum is fixed by the width of the observed CCF. To maintain orthogonality between jump parameters we replaced s with the full width at half maximum (FWHM) of the CCF, which is approximated by Here vg is the FWHM of the gaussian representing the local stellar and instrumental profile: The Markov chain is constructed by repeatedly computing χ 2 for a sequence of parameter sets {v sin I, λ, v ccf , γ}. At each step k, a set of trial parameters is computed by applying gaussian perturbations to the preceding set: v sin I k = v sin I k−1 + f σv sin I G(0, 1) The scale factor f is of order unity. The goodness-of-fit is recalculated for the new parameter set, and the change ∆χ 2 = χ 2 k − χ 2 k−1 is computed. The new set is either accepted or rejected according to the Metropolis-Hastings rule: if ∆χ 2 < 0 the set is accepted, and if ∆χ 2 > 0, the new step is accepted with probability e −∆χ 2 /2 . If the proposed parameters are accepted, they are written to the next record in the Markov chain. If not, the previous set of values is recorded. The scale factor f is tuned to ensure that roughly 25 percent of all steps are accepted once the chain has settled down to a steady state; by trial-and-error we found that f = 0.7 achieved the desired acceptance rate.
The chain evolved to a steady state within about 100 accepted proposals. We terminated the initial "burn-in" phase after at least 100 proposals had been accepted, when the value of χ 2 k exceeded the median χ 2 of all previous steps for the first time (Knutson et al. 2008). At this stage we refined our estimates of the variances of the binned CCF data, replacing the initial estimates σ 2 ij with σ 2 ij χ 2 min,j /n d . Here χ 2 min,j is the contribution of the jth spectrum to the χ 2 statistic for the best-fitting model in the chain so far, and  n d is the number of degrees of freedom associated with the spectrum. This ensures that χ 2 is scaled correctly, and that the posterior probability distributions of the fitted parameters have the correct variances.
We then ran the chain for a further 100 proposal acceptances, and re-evaluated the variances of the four fitting parameters from the chains themselves. After a further short burn-in period, we performed a production run of 10 4 accepted proposals. The correlation lengths of all four chains were computed as described by Tegmark et al. (2004), and found to be of order 10 steps in all four parameters.
The profile decomposition method carries higher computational overheads than the anomalous RV approaches developed by Winn et al. (2005) and Hirano et al. (2009). Substantial speed gains can be realised by performing the numerical integration to generate the synthetic stellar profile only once per step in the Markov chain, and shifting the profile to the appropriate orbital velocity at the time of each observation. The MCMC analyses presented here, with chains comprising 10 4 accepted steps, incurred execution times of order 50 minutes on a standard 2.66 GHz laptop processor, as opposed to a few minutes for a conventional analysis.

RESULTS
In Figs. 4 and 5 we show the correlation diagrams for all four jump parameters, for the nights of 2006 September 7/8 and 2007 August 28/29. In these plots we show the FWHM vg of the intrinsic gaussian profile, rather than the FWHM v ccf of the convolved profile. Contrary to our initial expectations, we find that v sin I and vg are only very weakly correlated. This confirms that the two independent measures of both quantities (from the shape of the stellar profile and the width and trajectory of the missing-light signature) serve to fix their values unambiguously.
The four fitted parameters for each of the two transits are listed in Table 2, for three different values of the linear limb-darkening coefficient u. These values bracket the value u 0.725 from the atlas9 models of Claret (2000) for a stellar effective temperature of 5050K, log g = 4.5 and nearsolar metallicity, as reported by Bouchy et al. (2005) for HD 189733. The tabulated values illustrate that the best-fitting values of v sin I and vg are only very weakly dependent on the precise value of u.
Although the two sets of observations were made nearly a year apart, the fitted values of the four parameters agree within their uncertainty ranges. The radial velocity γ of the system's centre of mass remained constant to within 10 m s −1 . The spin-orbit misalignment angle λ is found to be slightly negative in both years. We find an absolute value slightly, but not significantly, smaller than the (0.85 ± 0.30) • measured by Triaud et al. The uncertainty in λ is also 0.3 degrees, comparable with the uncertainty derived by Triaud et al. (2009) from the same data.
Our value of v sin I = 3.10 ± 0.03 km s −1 is significantly lower than the 3.32 +0.02 −0.07 km s −1 obtained by Triaud et al. (2009) using an uncorrected fit to the HARPS radial velocities. As Triaud et al. noted, the radial-velocity measurements obtained during transit by gaussian fitting show systematic departures from the intensity-weighted average velocity of the unobscured parts of the stellar disc. They modelled the departure of the radial-velocity centroid from the gaussian-fitting velocity as a function of v sin I, and found that v sin I = 3.05 km s −1 yielded the best fit to the observed pattern of residuals. An identical pattern of residuals was seen in the earlier R-M study of this system by Winn et al. (2006). Their empirical correction procedure yielded a lower v sin I = 2.97 ± 0.22 km s −1 . Our direct tracking of the missing-starlight component in the CCF eliminates the need for such empirical corrections to indirect measures of the missing-light velocity, and yields a more precise result that lies within the error ranges of both studies.

DISCUSSION AND CONCLUSIONS
In principle we could determine a lower limit on the stellar radius from v sin I and the photometric rotation period, which was determined by Henry & Winn (2008) to be 11.953 ± 0.009 days. The very close alignment of the projected orbital and spin axes in the plane of the sky makes it highly probable that the inclination I of the stellar spin axis to the line of sight should be very close to the inclination i = 85.5 degrees of the planet's orbit. Assuming sin I = sin i = 0.997, our value for v sin I yields a direct estimate of the stellar radius R * = 0.732 ± 0.007 R , independently of the usual assumptions involving isochrone fits to the stellar density and effective temperature. This is significantly smaller than the R * = 0.77 ± 0.01 R determined by Triaud et al. (2009), and marginally less than the R * = 0.75 ± 0.01 R found by Pont et al. (2007). We suspect that the discrepancy arises from differential rotation of the stellar photosphere. Gaudi & Winn (2007) pointed out that under some circumstances, the RM effect could be used to measure stellar differential rotation. As with most previous studies of this kind, we have measured v sin I under the assumption that the star rotates as a solid body. In all stars for which differential rotation has been measured directly, however, a solar-like pattern is found: the equator rotates faster than the poles, with the departure from the equatorial spin rate having a sine-squared dependence on latitude. The typical difference in rotational angular frequency ∆Ω = Ωequator − Ω pole is found to be of order 0.04 radian per day for rapidly-rotating stars K dwarfs similar to HD 189733, and only weakly dependent on rotation rate (Barnes et al. 2005). This is very similar to the solar value of ∆Ω.
If the starspots giving rise to the photometric modulation signal are concentrated at the same stellar latitude as the planet's path across the stellar disc, we should derive the correct value of the stellar radius. Unfortunately we have no way of determining the dominant latitude from which the modulation signal originates. If the orbital and spin axes are closely aligned, the impact parameter of the planet's path across the star carries it across a stellar latitude of 39 • in the hemisphere facing away from the observer. Although Pont et al. (2007) noted the passage of the planet across one or two isolated dark spots, these spots are significantly foreshortened. It is reasonable to expect the modulation signal to arise mainly from spots at intermediate to low latitude in the observer's hemisphere. A difference in latitude of 20 or 30 degrees is sufficient to give a difference in angular velocity δΩ 0.015 radians per day between the main spot belts and the latitudes traversed by the planet. This discrepancy yields a stellar radius that is too small by about 3 percent, which is sufficient to reconcile our radius estimate with those of Pont et al. (2007) and Triaud et al. (2009).
We conclude that direct modelling of the crosscorrelation function during a planetary transit is feasible even for a host star such as HD 189733 whose v sin I is comparable to the intrinsic line width. The profile decomposition method described here yields values and error estimates for the stellar spin rate and orbital obliquity that agree closely with the results of previous studies of this system. The parameter values are free of the systematic errors that occur when velocity measurements are derived from gaussian fits to line profiles with inherent time-variable asymmetry. This circumvents the need for semi-empirical corrections when modelling the Rossiter-McLaughlin effect. The method makes the fullest use of all the information present in the cross-correlation function: v sin I is tightly constrained by the shape of the stellar rotation profile as well as by the spectral signature of the starlight blocked by the planet. The stellar radius obtained from the photometric rotation period and the stellar v sin I agrees well with values obtained from high-precision transit photometry, if a modest amount of differential rotation is present.
The method described here could in principle be extended to observations for which an iodine cell is used to track the spectrograph PSF and wavelength scale. The analysis of such observations requires a reference spectrum of the host star, taken with the same instrument without the iodine cell. An artificial reference spectrum could be generated by convolving a limb-darkened rotation profile with the the spectrum of a narrow-lined star of the same spectral type. The Doppler signature of the light blocked by the planet could be mimicked by scaling, shifting and subtracting the same narrow-lined spectrum from the broad-ened stellar spectrum; the result would then be used as the reference spectrum. A very similar procedure was described by Winn et al. (2005) for calibrating the departure of the velocities measured during a transit of HD 209458b from the predictions of the Ohta, Taruya, & Suto (2005) model, using high-resolution solar spectra.