Improved Modeling of the Rossiter-McLaughlin Effect for Transiting Exoplanets
Teruyuki Hirano, Yasushi Suto, Joshua N. Winn, Atsushi Taruya, Norio Narita, Simon Albrecht, Bun'ei Sato
aa r X i v : . [ a s t r o - ph . E P ] A ug Improved Modeling of the Rossiter-McLaughlin Effect for TransitingExoplanets
Teruyuki Hirano , , Yasushi Suto , , , Joshua N. Winn , Atsushi Taruya , , , Norio Narita ,Simon Albrecht , and Bun’ei Sato [email protected] ABSTRACT
We present an improved formula for the anomalous radial velocity of the star dur-ing planetary transits due to the Rossiter-McLaughlin (RM) effect. The improvementcomes from a more realistic description of the stellar absorption line profiles, takinginto account stellar rotation, macroturbulence, thermal broadening, pressure broaden-ing, and instrumental broadening. Although the formula is derived for the case in whichradial velocities are measured by cross-correlation, we show through numerical simu-lations that the formula accurately describes the cases where the radial velocities aremeasured with the iodine absorption-cell technique. The formula relies on prior knowl-edge of the parameters describing macroturbulence, instrumental broadening and otherbroadening mechanisms, but even 30% errors in those parameters do not significantlychange the results in typical circumstances. We show that the new analytic formulaagrees with previous ones that had been computed on a case-by-case basis via numericalsimulations. Finally, as one application of the new formula, we reassess the impact ofthe differential rotation on the RM velocity anomaly. We show that differential rotationof a rapidly rotating star may have a significant impact on future RM observations.
Subject headings: planets and satellites: general – planets and satellites: formation –stars: rotation – techniques: radial velocities – techniques: spectroscopic Department of Physics, and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute ofTechnology, Cambridge, MA 02139 Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan Research Center for the Early Universe, School of Science, The University of Tokyo, Tokyo 113-0033, Japan Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544 Institute for the Physics and Mathematics of the Universe (IPMU), The University of Tokyo, Chiba 277-8582,Japan National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo, 181-8588, Japan Department of Earth and Planetary Sciences, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku,Tokyo, 152-8551, Japan
1. Introduction
Transiting exoplanetary systems provide valuable opportunities to learn about the natureof the exoplanets, their orbits, and their host stars. In particular, when we measure the radialvelocity (RV) of a star during a planetary transit, we see an anomalous Doppler shift (beyond theusual orbital RV) which is called the Rossiter-McLaughlin (hereafter, RM) effect (Rossiter 1924;McLaughlin 1924; Hosokawa 1953; Albrecht et al. 2007). It arises because a portion of the rotatingstellar disk is blocked by the planet. The partial occultation brings about a distortion in thespectral lines, which is manifested as an anomalous RV depending on the position of the planet onthe stellar disk (see, e.g., Queloz et al. 2000; Ohta et al. 2005; Winn et al. 2005; Narita et al. 2007;Collier Cameron et al. 2010). The time variation of the RV anomaly reveals the (sky-projected)angle λ between the planetary orbital axis and the stellar spin axis. Measurements of this anglehave proved to be an important observational clue to the origin of close-in giant exoplanets.It is widely assumed that close-in gas giants, of which more than 100 are known, formed at a fewAU away from their host stars and subsequently “migrated” inward (Lin et al. 1996; Lubow & Ida2010). Many planetary migration scenarios have been investigated, and some of them predict smallvalues of λ ≈ ◦ while others allow larger spin-orbit misalignment angles (e.g., Fabrycky & Tremaine2007; Wu et al. 2007; Nagasawa et al. 2008; Chatterjee et al. 2008). The observed distribution of λ and its dependence on the host star properties (such as masses and ages) can be important cluesto understand the origin of close-in giant exoplanets (Winn et al. 2010a; Fabrycky & Winn 2009;Morton & Johnson 2011).Observations of the RM effect have now become almost routine (Winn 2010; Moutou et al.2011). However, it is important to remember that the relationship between the observed RVanomaly, and the position of the planet on the stellar disk, is not completely straightforward. Thisis because the RM effect is actually a spectral distortion, even though it is frequently studiedas though it were a pure Doppler shift. Many alternatives have been pursued to calibrate therelationship between the observed signal and the underlying parameters of the planet and star.Ohta et al. (2005) and Gim´enez (2006) derived analytic formulas for the RV anomaly based on thecomputation of the first moment (the intensity-weighted mean wavelength) of distorted spectrallines. This approach is simple and convenient because the computed velocity anomaly does notdepend on the intrinsic shape of spectral lines (Hirano et al. 2010), and has been useful for quickcomputations where high accuracy is not essential, and for gaining insight into the parameterdependence of the RM velocity anomaly. Winn et al. (2005), however, noted that the analyticformula by Ohta et al. (2005) (hereafter, the OTS formula) deviates from the results based ona more realistic numerical calibration; they simulated spectra exhibiting the RM effect for manydifferent positions and sizes of the planet, and then analyzed the mock spectra with the samedata analysis codes that are used routinely to derive precise RVs with the High Resolution Echelle An alternative is to model the line profiles directly, as has been done by Albrecht et al. (2007) andCollier Cameron et al. (2010), which can be advantageous in some circumstances. §
5) is devoted to discussion and summary. 4 –
2. Derivation of the New Analytic Formula for the RM Effect
In this section, we derive the new analytic formula that describes the velocity anomaly duringa transit. We begin with our description of the stellar absorption lines. We follow the formulationby Hirano et al. (2010) but slightly change the basic equations in order to describe the stellar lineprofiles more realistically. Since the velocity field on the stellar surface is of primary importance, itis more convenient to express all the functions in terms of velocity rather than wavelength. In whatfollows, the velocity component v indicates the velocity shift relative to the center of an absorptionline. This is related to the wavelength shift ∆ λ by the usual formula ∆ λ/λ = v/c , where λ isthe central wavelength of the absorption line and c is the speed of light. Following the model ofspectral lines by Gray (2005), we write a stellar line shape F star ( v ) as F star ( v ) = − S ( v ) ∗ M ( v ) , (1)where S ( v ) is the intrinsic stellar line shape in the absence of stellar rotation and macroturbulence(for which we will give an explicit expression later), and M ( v ) is the broadening kernel due tostellar rotation and macroturbulence . The symbol ∗ indicates a convolution between two functions.Since the continuum level and the normalization factor in the spectrum do not affect the result inestimating the velocity anomaly during a transit, for convenience we subtract the continuum levelso that F star ( v ) becomes zero in the limit of v → ±∞ . Furthermore we normalize the spectrum sothat Z ∞−∞ F star ( v ) dv = − . (2)The minus sign indicates that F star ( v ) describes an absorption line. The rotational-macroturbulencebroadening kernel M ( v ) is calculated by disk-integrating the Doppler-shift component of the stellarsurface due to both stellar rotation and macroturbulence. We adopt “the radial-tangential model”for macroturbulence, for which the kernel in the absence of rotation isΘ( v ) = 12 √ π (cid:20) ζ cos θ e − (cid:16) vζ cos θ (cid:17) + 1 ζ sin θ e − (cid:16) vζ sin θ (cid:17) (cid:21) , (3)where ζ is the macroturbulent velocity parameter and θ is the angle between our line-of-sight andthe normal vector to the local stellar surface (Gray 2005, page 433). The angle θ is related to thecoordinate ( x , y ) on the stellar disk bycos θ = s − x + y R s , sin θ = p x + y R s , (4)where the y -axis is taken to be along the sky projection of the stellar spin axis, and R s is the radiusof the star. Assuming a quadratic limb-darkening law, the disk-integrated line broadening function We here assume a symmetric line profile and ignore the convective blueshift (CB) effect, discussed byShporer & Brown (2011). M ( v ) = Z Z entire disk − u (1 − cos θ ) − u (1 − cos θ ) π (1 − u / − u /
6) Θ( v − x Ω sin i s ) dx dyR s , (5)where u and u are the limb-darkening coefficients, Ω is the angular spin velocity of the star, and i s is the inclination angle of the stellar spin axis relative to the line of sight (Gray 2005). The Dopplershift − x Ω sin i s in the function Θ( v ) is the consequence of stellar rotation, neglecting differentialrotation. As Gray (2005) pointed out, the broadening kernel M ( v ) cannot be expressed as aconvolution of the two different broadening kernels of the stellar rotation and the macroturbulence.As we will show, the coupling between rotational broadening and macroturbulent broadening playsan important role in estimating the velocity anomaly due to the RM effect, especially when themacroturbulent velocity is appreciable when compared to the rotational velocity of the star (see thedifference between line profiles with and without macroturbulence shown in Figure 1). Indeed, thiscoupling between rotation and macroturbulence was neglected in the previous numerical calibrationsby Winn et al. (2005) and others. 6 – −20 −10 0 10 20velocity [km s −1 ]0.00.20.40.60.81.0 i n t en s i t y v sini = 0 km s −1 , ζ = 0 km s −1 v sini = 6 km s −1 , ζ = 0 km s −1 v sini = 6 km s −1 , ζ = 4 km s −1 intrinsicstellar rotationrotation + macroturbulence Fig. 1.— A schematic plot of the line profile during a planetary transit. Each line has a differentbroadening kernel. For visual clarity, the line profiles are vertically shifted by 0 . β = 1 km s − . Red line: after convolvingwith a pure-rotational broadening kernel (no macroturbulence), with v sin i s = 6 km s − . Blueline: after convolving with a rotational-macroturbulent broadening kernel with v sin i s = 6 km s − and ζ = 4 km s − . For the latter two cases (red and blue lines) the spectral contribution fromthe portion occulted by the planet has been subtracted from the profiles, assuming a planet with( R p /R s ) = 0 .
01. The line profile with macroturbulence (blue) has elongated wings and the transitsignal is nearly invisible.Next, we compute the line shape during a planetary transit. During a transit, the spectralcontribution of the portion blocked by the planet is written as F planet ( v ) = − S ( v ) ∗ M ′ ( v ) , (6)where M ′ ( v ) indicates a kernel similar to that given in Equation (5) but for which the disk-integration should only be performed over the blocked part of the stellar surface, rather than theentire stellar disk. As long as the planet is sufficiently small relative to the star, the Doppler shift − x Ω sin i s in Equation (5) is nearly constant over the integration region. Thus, if we define X asthe x -coordinate of the intensity-weighted center of the eclipsed portion of the star, we can remove 7 –the macroturbulence kernel Θ( v ) from the integral and define the following two useful quantities: f ≡ Z Z occulted portion − u (1 − cos θ ) − u (1 − cos θ ) π (1 − u / − u / dx dyR s , (7) v p ≡ X Ω sin i s , (8)so that F planet ( v ) becomes F planet ( v ) = − f S ( v ) ∗ Θ( v − v p ) . (9)The first quantity, f , is the instantaneous fractional decrease in flux due to the transit. The secondquantity, v p , is the rotational radial velocity of the occulted portion of the stellar disk, which isoccasionally referred to as the “subplanet velocity.” With these definitions the stellar line profileduring a transit is expressed as F transit ( v ) ≡ F star ( v ) − F planet ( v ) = − S ( v ) ∗ M ( v ) + f S ( v ) ∗ Θ( v − v p ) . (10)It should be noted that the macroturbulent kernel Θ( v ) remains in the modeled transit line profile.This treatment is necessary since the two effects of rotational broadening and macroturbulence arecoupled with each other. In short, the line profile during a transit expressed by Equation (10) isdifferent from the line profile modeled by Hirano et al. (2010) (Eq.[11]) in two senses: Equation(10) explicitly involves the effect of macroturbulence, and it is expressed in terms of velocity. 8 –Table 1: Summary of symbols used in this paper.Symbol Meaning Typical Range f the instantaneous fractional decrease in flux due to the transit (Eq. [7]) v p the subplanet velocity (Eq. [8]) ± v sin i s ∆ v the velocity anomaly due to the RM effect - u , u the limb-darkening parameters for the quadratic limb-darkening law v sin i s the stellar spin velocity - T eff the stellar effective temperature - i s the inclination of the stellar spin axis measured from our line-of-sight ◦ - 90 ◦ l the latitude on the stellar surface ± ◦ R s the stellar radius - x, y the position of the transiting planet on the stellar disk ± R s α the coefficient of differendtial rotation ± . β the Gaussian dispersion of spectral lines (Eq. [20]) 2.5 - 4.5 km s − γ the Lorentzian dispersion of spectral lines 0.5 - 1.5 km s − ζ the macroturbulence dispersion 2.0 - 6.5 km s − θ the angle between the line-of-sight and the normal vector to the stellar surface ◦ - 90 ◦ λ the spin-orbit misalignment angle ± ◦ ξ the microturbulence dispersion 0.0 - 2.0 km s − σ the frequency in the Fourier domain -Ω the angular velocity of the stellar spin -Armed with the preceding results, we now express the velocity anomaly ∆ v during a transitin terms of the fractional flux decrease f and the subplanet velocity v p . Basically, we followHirano et al. (2010) in order to compute the best-fit value for the anomalous RVs; they cross-correlated the spectrum during a transit with a stellar template spectrum, and then calculated thebest-fit value for the velocity anomaly ∆ v by maximizing the cross-correlation function C ( x ): dC ( x ) dx (cid:12)(cid:12)(cid:12) x =∆ v = 0 , (11) C ( x ) ≡ Z ∞−∞ F star ( v − x ) F transit ( v ) dv. (12)To proceed further, we need a specific model for the intrinsic line shape S ( v ). We here adopt theVoigt function for S ( v ): S ( v ) = V ( v ; β, γ ) ≡ G ( v ; β ) ∗ L ( v ; γ ) , (13) G ( v ; β ) ≡ β √ π e − v /β , (14) L ( v ; γ ) ≡ π γv + γ , (15) 9 –where β is the thermal velocity parameter and γ is the Lorentzian velocity parameter (due topressure broadening or natural broadening). These parameters are related to individual stellarproperties such as the effective temperature, surface gravity, and the nature of each absorptionline. Some line profiles of especially strong absorption lines (such as the Na D lines) are saturatedand intrinsically different from the Voigt function in shape. However, most of the lines in thewavelength region used in RV analyses are relatively weak, by design, and are well approximatedby the Voigt function in the absence of the stellar rotation and macroturbulence.Substituting Equations (1) and (10) into Equations (11) and (12), we compute the velocityanomaly ∆ v due to the RM effect. Since further calculations are mathematically complicated, wedescribe the detail of the derivation in Appendix A and write down the result alone:∆ v ≈ − f π Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) ˜Θ( σ ) sin(2 πσv p ) σdσ Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) n ˜ M ( σ ) − f ˜Θ( σ ) cos(2 πσv p ) o σ dσ , (16)where ˜ M ( σ ) ≡ Z ∞∞ M ( v ) e − πiσv dv = Z − u (1 − √ − t ) − u (1 − √ − t ) − u / − u / × { e − π ζ σ (1 − t ) + e − π ζ σ t } J (2 πσv sin i s t ) tdt, (17)˜Θ( σ ) = 12 (cid:2) exp {− ( πζ cos θ ) σ } + exp {− ( πζ sin θ ) σ } (cid:3) , (18)where J n ( x ) is the Bessel function of the first kind (see also Appendix B). Equation (16) is themain finding in the present paper and will be used in the comparison with simulated results.
3. Validity of the Analytic Formula3.1. Comparison with Numerical Simulations for Subaru/HDS
One of the major differences between the derivation of the analytic formula (Eq. [16]) and themanner in which RV data are actually analyzed is that the analytic formula is based on the cross-correlation method while the data analysis (for Subaru/HDS and Keck/HIRES at least) adoptsthe forward modeling method using the Iodine cell for a precise wavelength calibration. Anotherdifference between them is that the analytic formula assumes only one absorption line in its deriva-tion while the actual stellar spectra have many lines differing in depth, shape, and so on. Thus, inorder to test the validity of Equation (16), we perform the mock data simulation described belowand compare the results with the analytic formula. 10 –
The mock data simulation with the Iodine RV calibration is described in detail by Winn et al.(2005) and Narita et al. (2009a). In order to generate mock spectra during a transit, we begin withtheoretically synthesized spectra by Coelho et al. (2005), whose spectral lines are modeled by in-corporating thermal broadening (including microturbulence), and Lorentzian (natural or pressure)broadening. Since intrinsic line profiles depend on the stellar type, we first test the comparison fora G0 star, whose effective temperature is 6000 K. When we obtain the synthetic spectrum from thesynthetic spectrum library, we assume a representative value of the surface gravity (log g = 4 .
0) fora G0 star with a planet and the solar abundance for metallicity ([Fe/H]=0.0).Using the synthetic spectrum, we generate the mock spectra during a transit by the followingprocedure.1. We broaden the synthetic spectrum by convolving the rotational-macroturbulent kernel M ( v )(Eq. [5]), so that it represents the template spectrum of the actual star. We adopt thelimb-darkening parameters u = 0 .
43 and u = 0 .
31 (Claret 2004), and the macroturbulenceparameter ζ = 4 . − (Valenti et al. 2005). We try three different values for the rotationalvelocity of the star: v sin i s = 2 . − , 5.0 km s − , and 7.5 km s − .2. We Doppler-shift the original unbroadened spectrum by v p , multiply by f , and then convolvethe spectrum with the macroturbulence kernel Θ( v ) (Eq. [3]), so that the resultant spectrumrepresents the spectral contribution from the portion occulted by the transiting planet. Foreach of the three values of v sin i s , we consider f = 0 . v p evenly spaced from − v sin i s to + v sin i s , yieldingin total 105 different points in the ( f, v p ) grid for each value of v sin i s . The macroturbulentbroadening depends not only on the x -coordinate but also on the y -coordinate on the stellardisk (see Eqs. [3] and [4]). For simplicity, however, we assume y = 0 when we generate themock transit spectra.3. We generate the mock transit spectra by subtracting the Doppler-shifted and scaled spectrain the second step from the broadened spectra created in the first step.4. We multiply the mock transit spectra by the iodine transmission spectrum used for cali-bration, and convolve the Star+I spectrum with the representative instrumental profile ofSubaru/HDS for the case of the slit width being 0.8 ′′ , corresponding to the spectral resolutionof R ∼ .We then take these 105 mock spectra for each of three values of v sin i s , and use them as inputsto the RV analysis routine for Subaru/HDS. The RV analysis using the iodine cell is described in The choices for the instrumental profiles used in the mock simulation are explained later.
11 –detail by Sato et al. (2002) (for Subaru/HDS) and Butler et al. (1996) (for Keck/HIRES). For eachpoint of the ( f, v p ) grid, the RV analysis routine outputs a velocity anomaly ∆ v due to the RMeffect. We compare the results ( f, v p , ∆ v ) based on the mock data simulation with the analytic ex-pression (16). Figure 2 shows the comparison between the simulated data points ( f, v p , ∆ v ) andthe analytic formula (Eqn. [16]). The three different panels are for the three different values of v sin i s : (top) 2.5 km s − , (middle) 5.0 km s − , and (bottom) 7.5 km s − . For each of the fivevalues of f , the data points indicated by color symbols show the simulated results; f = 0 .
004 inblack, f = 0 .
008 in red, f = 0 .
012 in blue, f = 0 .
016 in purple, and f = 0 .
020 in green.In order to draw the analytic curves, we need to choose values of β and γ . Although they arerelated to intrinsic stellar line profiles, each spectral line has its own values for β and γ , so we needto know the “effective” values of those line parameters in order to make a comparison betweenthe analytic formula and the simulated results. For the purpose, we employ an autocorrelationmethod. By autocorrelating the synthetic spectrum for a G0-type star, we obtain an effective lineprofile of the spectrum (see Appendix C for details). Since there is a strong degeneracy between thetwo stellar line parameters β and γ , we fix β based on a simple physical principle and estimate theLorentzian dispersion γ from the effective line profile. In principle, the intrinsic Gaussian dispersion β in each spectral line is determined by the effective temperature T eff of the star and the “themicro-turbulence” dispersion ξ as β = s k B T eff µ + ξ , (19)where k B and µ are the Boltzmann constant and the mass of the atom (or molecule) in question,respectively (Gray 2005). In addition, as we will show later, the simulated velocity anomalies alsodepend on the instrumental profile which we assume in generating mock spectra. Therefore, wehere adopt an ad hoc assumption that the Gaussian width parameter β in Equation (16) dependsalso on the width of the instrumental profile so that β is expressed as β = q β + β = s k B T eff µ + ξ + β , (20)where β IP is the dispersion of Gaussian which best-fits the representative instrumental profileadopted in observations.Substituting T eff = 6000 K = 0 .
517 eV, µ = 52 GeV/c (the mass of an iron atom ), ξ = 1 . − (Coelho et al. 2005), and β IP = 3 . − (the dispersion of the instrumental profile Iron lines are most numerous in the wavelength range used for the RV analysis with the iodine cell.
12 –of Subaru/HDS we assumed in generating mock transit spectra, corresponding to the spectralresolution of R ∼ β = 1 . − and β = 4 . − . Fitting the effective lineprofile with the Voigt function assuming β = 1 . − , we obtain γ = 0 . − (AppendixC). In each panel (different v sin i s ) of Figure 2, the analytic curves based on Equation (16) areshown in the same colors as the simulated results indicated by symbols for each value of f . Weadopt ( β, γ )=(4.0 km s − , 0.9 km s − ). The other parameters ( u , u , v sin i s , ζ ) in Equation (16)are fixed at the same values used to make mock transit spectra. 13 – -60-40-20 0 20 40 60 -3 -2 -1 0 1 2 3 ∆ v [ m s - ] v p [km s -1 ]G0-type Starv sin i s = 2.5 km s -1 f = 0.004f = 0.008f = 0.012f = 0.016f = 0.020-100-80-60-40-20 0 20 40 60 80 100 -6 -4 -2 0 2 4 6 ∆ v [ m s - ] v p [km s -1 ]G0-type Starv sin i s = 5.0 km s -1 f = 0.004f = 0.008f = 0.012f = 0.016f = 0.020-150-100-50 0 50 100 150 -8 -6 -4 -2 0 2 4 6 8 ∆ v [ m s - ] v p [km s -1 ]G0-type Starv sin i s = 7.5 km s -1 f = 0.004f = 0.008f = 0.012f = 0.016f = 0.020 Fig. 2.— Simulated velocity anomalies due to the RM effect v.s. the analytic formula (Eq. [16]) asa function of the subplanet velocity v p in the case of a G0-type star: v sin i s = 2 . − (top), 5.0km s − (middle), and 7.5 km s − (bottom). In each panel, symbols in black, red, blue, purple, andgreen indicate the simulated data points for f = 0 . v sin i s is small, the velocity anomaly curves are nearly linearand can be approximated as ∆ v ≈ − f v p . As Winn et al. (2005) and other authors have pointedout, however, they are significantly curved for larger values of v sin i s . In the forward modeling method using the iodine cell, radial velocities v RV of a star arecomputed by modeling each observed spectrum I obs ( λ ) with the following equation: I obs ( λ ) = k [ A ( λ ) T ( λ (1 − v RV /c ))] ∗ IP , (21)where A ( λ ) and T ( λ ) are the transmission spectrum of the Iodine cell and the template spectrum ofthe same star, respectively (Sato et al. 2002; Butler et al. 1996). Outside of a transit, the intrinsicstellar line profile of each spectrum is supposed to be the same as that of the template spectrum,and therefore v RV is not affected by the instrumental profile which may be variable during anobservation. During a transit, however, the intrinsic stellar line profile is distorted due to thepartial occultation, and it is not clear if the IP affects the radial velocity v RV including the RMvelocity anomaly ∆ v . Thus, in order to see if the instrumental profile used in the mock data analysisaffects the RM velocity anomaly, we repeated the same mock data analysis described above, butwith a different instrumental profile. In the simulation above, we have fixed the instrumental profileso that it corresponds to the spectral resolution of R = 45000 for Subaru/HDS. This time, we adoptthe instrumental profile with the spectral resolution of R = 90000. In this case, the dispersion ofGaussian fits the instrumental profile is approximately β IP = 1 . − . After generating manymock spectra for v sin i s = 2 .
5, 5.0, and 7.5 km s − as in Section 3.1.1, we analyzed them with theusual RV routine to obtain the velocity anomaly ∆ v .In order to quantify the differences between the two cases of the instrument profile, we computethe following statistics for each case of the instrumental profiles: D ( β, γ ) ≡ vuut X f,v p ,v sin i s (cid:12)(cid:12)(cid:12) ∆ v sim ( f, v p ) − ∆ v ana ( f, v p , β, γ ) f · v sin i s (cid:12)(cid:12)(cid:12) , (22)where ∆ v sim ( f, v p ) is the simulated velocity anomaly for each point of the ( f, v p ) grid, and ∆ v ana ( f, v p , β, γ )is the value computed by Equation (16) as a function of f , v p , β , and γ . When we compute∆ v ana ( f, v p , β, γ ), all the parameters other than β and γ are fixed at exactly the same values as areadopted in the numerical simulation (such as ζ and u , u ). The summation in the above expressionis performed over the 315 (= 105 ×
3) data points on the ( f, v p ) grid with three different values of v sin i s . The statistics, D ( β, γ ), indicates the degree of agreement between the analytic formula andthe simulated results in terms of the representative velocity anomaly f v sin i s . For each of the instru-mental profiles ( R = 45000 and R = 90000), we compute D ( β, γ ), for 2 . − ≤ β ≤ . − and 0 . − ≤ γ ≤ . − . 15 – γ [ k m s - ] β [km s -1 ]R =90000R = 45000 Fig. 3.— The contour of D ( β, γ ) for a G0-type star. The regions surrounded by the two blueand the two red curves are the best-fit regions of the analytic formula, where D ( β, γ ) ≤ .
005 for R = 45000 and R = 90000, respectively. The values of ( β, γ ) estimated by the line analysis areshown by the blue and red crosses for R = 45000 and R = 90000, respectively.As a result of computing D ( β, γ ) in the ( β, γ ) grid, we find the lowest value of D ( β, γ ) isapproximately 0.0045 for both of the two different instrumental profiles, representing very goodagreement. For instance, the dispersion of the simulated velocity anomalies around the analyticformula is expected to be less than ∼ .
25 m s − in case of f ≈ .
01 and v sin i s ≈ . − .This deviation is much smaller than the usual precision with which RVs can be measured.In Figure 3, we show contours of the goodness-of-fit statistic D ( β, γ ). The regions betweenthe two blue curves and two red curves are where D ( β, γ ) ≤ .
005 for R = 45000 and R = 90000,respectively. Since the Gaussian and Lorentzian dispersions β and γ strongly degenerate, thevalues of β and γ that fit the simulated results well are expected to be located in extended areas.Thus, we show the “best-fit regions” where the analytic formula agrees well with the simulatedresults. From Figure 3, it is obvious that the best-fit region is shifted toward smaller values of( β, γ ) when we adopt the higher spectral resolution. This result implies that the instrumentalprofile in spectroscopic observations does affect the velocity anomaly due to the RM effect. Weinterpret the results as follows; in the forward-modeling fitting procedure, the instrumental profile(a finite spectral resolution) plays a similar role as the other physical broadening mechanisms forspectral lines. For reference, we show ( β, γ )=(4.0 km s − , 0.9 km s − ) by the blue cross, whichare the intrinsic line parameters for R = 45000 estimated by the spectral line analysis in AppendixC. Also, substituting β IP = 1 . − into Equation (20) for the case of R ∼ β = 2 . − , which, along with γ = 0 . − (the same value as used in the comparison for R = 45000), is shown in Figure 3 by the red cross. It should be emphasized that the values of ( β, γ )shown by the blue and red crosses are estimated in a way independent of the mock data simulation(Appendix C) but are consistent with the regions where the analytic formula best agrees with thesimulated results, for the two different instrumental profiles.In summary, as a result of trying two different cases of the spectral resolution ( R ∼ β by adding theinstrumental broadening β IP in quadrature, which justifies the treatment shown in Equation (20). So far, we have considered a G0 star. In order to make sure that our analytic formula isapplicable for a variety of different stellar types, we consider stars with the effective temperaturesof T eff = 6500 K (an F5 star) and 5000 K (a K0 star) using the theoretically synthesized spectra,just as we did for a G0 star. F5 star
In generating mock spectra for an F-type star ( T eff = 6500 K), we broaden the syntheticspectrum of an F5 star (Coelho et al. 2005) assuming the rotational velocity of v sin i s = 5 . − . The other adopted parameters include u = 0 .
32 and u = 0 .
36 for thelimb darkening parameters (Claret 2004), and ζ = 6 . − for the macroturbulent velocityparameter (Gray 2005). Except for those parameters, we perform exactly the same simulationdescribed in Section 3.1.1. The simulated results for v sin i s = 10 km s − are shown in theupper panel of Figure 4, along with the analytic formula (Eq. [16]) assuming ( β, γ ) = (4 . − , 0.9 km s − ). These Gaussian and Lorentzian dispersions are also estimated by thecombination of Equation (20) adopting T eff = 6500 K, and the spectral line analysis usingthe auto-correlation method. Again, the analytic curves well describe the behavior of thesimulated results. They are also in good agreement with each other for v sin i s = 5 . − and 15 km s − , although for brevity we do not show all of those results here. 17 – -200-150-100-50 0 50 100 150 200 -10 -5 0 5 10 ∆ v [ m s - ] v p [km s -1 ]F5-type Starv sin i s = 10.0 km s -1 f = 0.004f = 0.008f = 0.012f = 0.016f = 0.020-60-40-20 0 20 40 60 -3 -2 -1 0 1 2 3 ∆ v [ m s - ] v p [km s -1 ]K0-type Starv sin i s = 3.0 km s -1 f = 0.004f = 0.008f = 0.012f = 0.016f = 0.020 Fig. 4.— The comparison between simulated results (symbols) and Equation (16) (solid curves)for an F5 star with v sin i s = 10 . − ( upper ) and for a K0 star with v sin i s = 3 . − ( lower ), respectively. K0 star
We also make the mock transit spectra for a K0 star ( T eff = 5000 K) and put theminto the RV routine. Since most of the K-type dwarfs are relatively slow rotators, we adoptsmall rotational velocities: v sin i s = 1 .
5, 3.0, and 4.5 km s − . We fix the limb-darkeningparameters and the macroturbulent velocity parameter to be u = 0 . u = 0 .
14 (Claret2004), and ζ = 2 . − (Valenti et al. 2005). We compare the simulated velocity anomalieswith Equation (16). As an example we show the result for v sin i s = 3 . − in the lowerpanel of Figure 4, for which we assume β = 3 . − and γ = 1 . − , as estimated byanalyzing the synthetic line profiles for G0 and F5 stars. 18 –All the numerical simulations show that our new analytic formula reproduces the simulated resultswithin the current RV precisions. The agreement also suggests that we can validate the previ-ously reported empirical relations for Subaru/HDS (i.e., Narita et al. 2009a,b, 2010a,b; Hirano et al.2011), which are based on the similar mock data simulations. In the previous subsection, we have shown that the analytic formula (Eq. [16]) gives an accuratedescription of the simulated velocity anomalies during a transit as long as we adopt appropriatevalues of β and γ for a given type of star. However, there are several practical issues that must beaddressed before the analytic formula is applied to real data analysis. First, the macroturbulentvelocity parameters that we assumed in both the application of the analytic formula and in themock data simulations will not be known a priori for real stars. Second, we have used theoreticallysynthesized spectra to generate mock transit spectra, but actual intrinsic line profiles in observedspectra may differ from theoretical ones. Finally, although we have adopted representative instru-mental profiles in both of the analytic formula and simulation, the instrumental profile is generallydependent on many factors such as temperature variations, the position on CCD, etc, even if weadopt the same spectrograph setups for observations. Therefore, it is important to investigate thesensitivity of Equation (16) to those parameters affecting the line profile (including the instrumen-tal profile). In this subsection, we perturb the values of those line parameters and examine theresulting changes to the outputs of the analytic formula. First we consider changes in the macroturbulence parameter. The macroturbulence dispersion ζ for a G0-type star like HD 209458 is empirically estimated as ζ = 4 . − (Valenti et al. 2005).We here change it by ±
30% with respect to that value (i.e. ζ = 3 . − ) and plot thevelocity anomaly curves for HD 209458 assuming the other line parameters as ( β, γ ) = (4 . − ,0.9 km s − ) and the stellar rotational velocity of v sin i s = 4 . − . The top panel in Figure5 shows the three cases of the macroturbulence dispersion: ζ = 4 . − (black), 3.0 km s − (red),and 5.6 km s − (blue). 19 – -4-2 0 2 4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-40-20 0 20 40 ∆ v [ m s - ] HD 209458v sin i s = 4.5 km s -1 ζ = 4.3 km s -1 ζ = 3.0 km s -1 ζ = 5.6 km s -1 -10-5 0 5 10 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-120-100-80-60-40-20 0 20 40 ∆ v [ m s - ] XO-3v sin i s = 18.5 km s -1 ζ = 6.2 km s -1 ζ = 4.3 km s -1 ζ = 8.1 km s -1 -4-2 0 2 4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-30-20-10 0 10 20 30 ∆ v [ m s - ] TrES-1v sin i s = 1.5 km s -1 ζ = 2.8 km s -1 ζ = 2.0 km s -1 ζ = 3.6 km s -1 Fig. 5.— Variations of RM velocity anomaly curves for HD 209458 (top), XO-3 (middle), and TrES-1 (bottom) based on Equation (16). In each panel, we change the macroturbulence parameter ζ by ±
30% (red and blue curves) from the literature-based values (black curves). The residuals ofthe red and blue curves from the black curves are also shown at the bottom for each system. Forthe spin-orbit misalignment angle λ , we assumed λ = 0 ◦ for HD 209458 and TrES-1, and λ = 37 . ◦ (Winn et al. 2009), respectively. 20 –Inspection of Figure 5 shows that the RV difference between the curves is at most ∼ − ( ∼ .
5% of the representative velocity anomaly ( R p /R s ) v sin i s for HD 209458), which is less thanthe RV precision from an observational point of view. Therefore, the choice of ζ for this type ofstar does not significantly affect the results.We also show the comparisons for other two systems: XO-3 and TrES-1. For an F5 star XO-3,we adopt β = 4 . − , γ = 0 . − , and v sin i s = 18 . − , and try three different casesof ζ : 4.3, 6.2, and 8.1 km s − (changed by ±
30% around the central value). Likewise, we computethe RM velocity anomaly for TrES-1, assuming β = 3 . − , γ = 1 . − and v sin i s = 1 . − . We choose values for the macroturbulence dispersion of ζ = 2 .
0, 2.8, and 3.6 km s − . Theresults are shown in the middle and bottom panels in Figure 5. The RV differences of the coloredcurves from the black ones are at most ≈ − ( ∼ .
5% of ( R p /R s ) v sin i s ) for XO-3 and ∼ − ( ∼ .
4% of ( R p /R s ) v sin i s ) for TrES-1, respectively. These values are less than the RVprecision that has been achieved for each system, which are approximately 8 m s − (Winn et al.2009) and 10 m s − (Narita et al. 2007), respectively. β and γ Next, we focus on the Gaussian and Lorentzian components of line profiles. Fixing the valuesof macroturbulence parameter ζ , we change β and γ in Equation (16), and plot the RM curvesfor the three systems above. Figure 6 presents the three different cases of ( β, γ ) for HD 209458.We change both of β and γ by ±
30% from the values estimated by the spectral line analysis(Appendix C) while we fix ζ at 4.3 km s − . According to Figure 6, the largest discrepancy of thetwo colored curves from the black one is . − ( ∼ .
1% of ( R p /R s ) v sin i s ). Although this isstill comparable to or less than the RV precision for this type of stars, it is slightly larger than thedifference in the curves for different ζ (see the top panel of Figure 5). 21 – -4-2 0 2 4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-40-20 0 20 40 ∆ v [ m s - ] HD 209458v sin i s = 4.5 km s -1 ( β , γ ) = (4.0, 0.9)( β , γ ) = (2.8, 0.6)( β , γ ) = (5.2, 1.2) Fig. 6.— The variation of RM velocity anomaly curve for HD 209458. We change the Gaussianand Lorentzian dispersion ( β, γ ) in Equation (16) by ±
30% (red and blue curve) from the valuesbased on the line analysis (black curve). The residuals of the red and blue curves from the blackone are also shown at the bottom.We also compare the analytic RM curves for different values of ( β, γ ) in case of XO-3 (F5-typestar) and TrES-1 (K0-type star). As a result, the differences in the RM curves are at most . − ( ∼ .
1% of ( R p /R s ) v sin i s ) for XO-3 and . . − for TrES-1 ( ∼ .
0% of ( R p /R s ) v sin i s ).These deviations are less than the current RV precision for each type of star. For a rapidly rotatingstar such as XO-3, the line profile is mainly determined by the stellar rotation and macroturbulence,and thus the intrinsic line profiles (thermal and Lorentzian broadening) are likely to be less impor-tant. On the other hand, for slowly rotating late-type stars as TrES-1, the intrinsic line parameters( β, γ ) are important for describing the line shapes. However, since the RM velocity anomaly is wellapproximated as ∆ v ≈ − f v p for narrow line profiles (slow rotators) and the velocity amplitudedue to the RM effect is comparably small, an inaccurate estimation for ( β, γ ) or a variation of theinstrumental profiles (less than ∼ v sin i s &
10 km s − ), themacroturbulence dispersion is more important than thermal and natural profiles, while the Gaussianand Lorentzian dispersions become more significant as for moderately rotating stars (3.0 km s − . v sin i s .
10 km s − ). Neither of the effects is important for slowly rotating stars ( v sin i s . . − ). 22 – In Section 3.1, we compared the analytic formula with simulated results based on the analysisroutine for Subaru/HDS. We have confirmed that Equation (16) well approximates the simulatedresults with deviations less than ∼ .
5% of the typical velocity anomaly scale (i.e. f v sin i s ). Itis very interesting to see if the new formula (Eq. [16]) also agrees with the simulated results inpreviously published papers using some instrument other than Subaru/HDS. If they are shown tobe consistent with each other, it suggests we no longer need to perform numerical case-by-casemock data simulations for RM analyses. If they do not agree, it may justify further investigationinto which approach is more accurate. Here, we consider four systems for which the RM effect hasbeen measured with Keck/HIRES: HAT-P-4, HAT-P-14, XO-3, and HAT-P-11. We compare Equa-tion (16) with the published empirical relations based on mock data simulations for Keck/HIRES.As discussed in Section 3.1.3, the instrumental profile affects the velocity anomaly due tothe RM effect. Therefore, in order to make a comparison between Equation (16) and the empiricalformulae, we need to know the representative dispersion of the instrumental profile of Keck/HIRES.We here adopt β IP = 2 . − in terms of velocity (a representative dispersion for the spectralresolution of R = 65000). Since other parameters are intrinsic stellar properties, we adopt similarvalues of β , γ , and ζ to the ones adopted for the comparison in Section 3.1 to draw the analyticRM curves.Table 2: Calibrations of the RM effect drawn from the literature, and choices of the stellar lineparameters used in Figure 7. The stellar parameters ( β, γ, ζ, v sin i s ) are expressed in km s − .System Empirical Relation β γ ζ v sin i s ReferenceHAT-P-4 ∆ v = − f v p [1 . − . (cid:16) v p v sin i s (cid:17) ] 3.2 0.9 4.1 5.5 Winn et al. (2011)HAT-P-14 ∆ v = − f v p [1 . − . (cid:16) v p v sin i s (cid:17) ] 3.3 0.9 6.5 8.4 Winn et al. (2011)XO-3 ∆ v = − f v p [1 . − . (cid:16) v p v sin i s (cid:17) ] 3.3 0.9 6.2 18.5 Winn et al. (2009)HAT-P-11 ∆ v = − f v p -10-5 0 5 10 -2 -1 0 1 2 R e s i dua l [ m s - ] time [hour]-40-20 0 20 40 ∆ v [ m s - ] (a)HAT-P-4v sin i s = 5.5 km s -1 λ = - 4.9 degEquation (16) ( ζ = 4.1 km s -1 )Keck SimulationGaussian formulaOTS formula -10-5 0 5 10-1.5 -1 -0.5 0 0.5 1 1.5 R e s i dua l [ m s - ] time [hour]-10 0 10 20 30 40 ∆ v [ m s - ] (b) HAT-P-14v sin i s = 8.4 km s -1 λ = 187.8 degEquation (16) ( ζ = 6.5 km s -1 )Keck SimulationGaussian formulaOTS formula -30-20-10 0 10 20 30 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-100-50 0 50 ∆ v [ m s - ] (c)XO-3v sin i s = 18.5 km s -1 λ = 37.3 degEquation (16) ( ζ = 6.2 km s -1 )Keck SimulationGaussian formulaOTS formulaHD3861 + rotation kernel -0.4-0.2 0 0.2 0.4 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-10-8-6-4-2 0 2 4 ∆ v [ m s - ] (d) HAT-P-11v sin i s = 1.5 km s -1 λ = 103 degEquation (16) ( ζ = 2.5 km s -1 )Keck SimulationGaussian formula Fig. 7.— The comparison between Equation (16) (black) and the empirical relations in literatures(red) for (a) HAT-P-4, (b) HAT-P-14, (c) XO-3, and (d) HAT-P-11. For reference, analytic curvesbased on the Gaussian formula by Hirano et al. (2010) and the OTS formula are plotted in blueand green, respectively. The residuals of the empirical curves and the previous analytic formulaefrom Equation (16) are also shown at the bottom in each panel.Figure 7 shows the comparison between Equation (16) and the empirical relations in literaturesfor the four systems. The empirical relations in the literatures are summarized in the second columnof Table 2 and shown with red curves in Figure 7. Analytic curves (Eq. [16]) are shown in black. 24 –The adopted stellar line parameters to draw the analytic curves are also summarized in Table 2.In Table 2, β , γ , ζ , and v sin i s are expressed in km s − . For reference, we plot the Gaussianformula by Hirano et al. (2010) (i.e. Eqn.[27]) in blue and the OTS formula in green, which isexpressed as ∆ v = − f v p / (1 − f ). In drawing the Gaussian formula, we employ the width β p ofthe intrinsic spectral line (in the absence of stellar rotation) approximated as a single Gaussian by β p = p β + ζ , where β and ζ are the values shown in Table 2.Regarding HAT-P-4, HAT-P-14, and HAT-P-11, for which the rotational velocities are rela-tively small, Equation (16) and the numerical calibration formulas (red) agree within a few m s − .For the HAT-P-11 system in particular, the two curves are almost indistinguishable. This is becauseEquation (16) is well approximated by ∆ v ≈ − f v p for the case in which stellar lines are sufficientlynarrow. On the other hand, the discrepancy between them for the case of XO-3 is more significant( ∼
10 m s − ) in comparison with other systems although they are in much better agreement incomparison to the previous analytic curves based on the Gaussian and the OTS formulae.There are at least two possible explanations for the discrepancies between Equation (16) andthe Keck simulated result. First, the line parameters ( β, γ, ζ ) adopted in Equation (16) to plotthe analytic curves in Figure 7 may be significantly different from the true values for the system.Specifically, the macroturbulence dispersion ζ is not well known for massive stars, and differentvalues of ζ lead to different results as we have shown in Section 3.2.1. In Section 3.1, on the otherhand, we have compared Equation (16) with simulated results adopting the same values of ζ inEquation (16) as the ones used to generate mock transit spectra.A second possible explanation for the discrepancy between the analytic formula and the Kecksimulation is that the coupling between the stellar rotation and macroturbulence (which was ignoredin the numerical calibrations) has a non-negligible impact on the results. In most cases, the mockdata simulations for Keck/HIRES used real observed spectra to create mock spectra during atransit. For example, to create mock spectra during a transit of the XO-3 system, Winn et al.(2009) began with the template spectrum of HD 3861, which is the same spectral type as XO-3 buthas a smaller rotational velocity (F5V, v sin i s = 2 . − ). Then, they broadened that spectrumwith a pure-rotational kernel to mimic the spectrum of XO-3. This is a good choice in the sensethat the starting spectrum already involves the macroturbulence, which is not well known for thosetypes of stars. However, as we have emphasized, the rotational and macroturbulence broadeningsare coupled and, strictly speaking, cannot be applied sequentially. To be more precise, we shouldfirst deconvolve the spectrum with the macroturbulence before convolving with a kernel includingthe effects of both rotation and macroturbulence.In order to test the second hypothesis, we perform another mock data simulation for XO-3.We begin with the F5-type synthetic spectrum and broaden it with M ( v ) in which we assume v sin i s = 2 . − and ζ = 6 . − , so that the resulting spectrum mimics the templatespectrum of HD 3861 (let us call it T HD3861 ( λ )). Then, we additionally broaden T HD3861 ( λ ) withthe pure-rotational broadening kernel without macroturbulence assuming v sin i s = 18 . − (we 25 –call the resulting spectrum as T XO3 , Keck ( λ )). This process is expected to reproduce the steps takenby Winn et al. (2009) in creating the mock spectrum of XO-3, although our starting spectrum is atheoretically synthesized one rather than that of HD 3861. By changing ( f, v p ), we generate manymock transit spectra from T XO3 , Keck ( λ ) and put the mock spectra into the RV analysis routine.As a result of fitting the output velocity anomaly ∆ v as a function of f and v p , we obtain thefollowing empirical relation:∆ v = − f v p h . − . (cid:16) v p . − (cid:17)i ± σ fit , σ fit = 0 . − , (23)where σ fit is the dispersion of the residuals of the simulated ∆ v from the best-fit curve. We plotthis empirical relation in Figure 7 (c) in purple. Although we still see differences between the Kecksimulation (red) and our new simulation (purple), they are reduced, and the new empirical relationfollows more closely the variation seen in the published numerical calibration. This suggests thatpart of the discrepancy (and perhaps most of the discrepancy) between Equation (16) and the Keckempirical formula is ascribed by the coupling between the stellar rotation and macroturbulencewhich was neglected in performing the numerical calibration. Note that the RV difference betweenthe two simulations (red and purple) is comparable to the dispersion of our simulated results( σ fit = 7 . − ).
4. Impact of Differential Rotations of Stars
One of the important applications of the new analytic formula is describing the stellar dif-ferential rotation in RM measurements. Gaudi & Winn (2007) previously pointed out the impactof differential rotations via RM measurements is negligible since they are beyond the precision ofthe RV measurements. However, their argument on the detectability is focused on relatively slowrotators. Since the RM effect is now measured for fairly rapid rotators ( v sin i s &
10 km s − ), it isimportant to see if the impact of stellar differential rotations is still negligible based on the precisemodeling of the RM effect described in this paper. In this section, we describe the impact of stellardifferential rotations in the velocity anomaly curves.Following Reiners (2003), we model the stellar angular velocity Ω as a function of the latitude l on the stellar surface ( l = 0 at the stellar equator) asΩ( l ) = Ω eq (1 − α sin l ) , (24)where Ω eq is the angular velocity at the equator. The coefficient α is the parameter describingthe degree of differential rotation and approximately 0.2 for the case of the Sun. If the transitingplanet is located at ( x, y ) on the stellar disk, its latitude l is estimated bysin l = yR s sin i s + s − x + y R s cos i s . (25) 26 –When the host star is differentially rotating, we must replace the constant Ω in Equations (8)and (B7) with Equation (24). In this case, the degeneracy between Ω and sin i s can be solved; forrigidly rotating systems, the velocity anomaly ∆ v depends solely on Ω sin i s , but with a differentialrotation, the specific choice of the stellar inclination i s also affects the result. -15-10-5 0 5 10 15 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-100-50 0 50 ∆ v [ m s - ] XO-3v sin i s = 18.5 km s -1 α = 0.0, i s = 90 deg α = 0.1, i s = 90 deg α = 0.2, i s = 90 deg α = 0.3, i s = 90 deg α = 0.4, i s = 90 deg -15-10-5 0 5 10 15 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 R e s i dua l [ m s - ] time [hour]-100-50 0 50 ∆ v [ m s - ] XO-3v sin i s = 18.5 km s -1 α = 0.0, i s = 90 deg α = 0.2, i s = 90 deg α = 0.2, i s = 60 deg α = 0.2, i s = 30 deg Fig. 8.— The RM velocity anomaly curves for the XO-3 system with and without differentialrotation. The blacks lines in each panel indicate the case when the host star is rigidly rotatingwhile the other colors are for the cases that the star is differentially rotating. In the left panel,the stellar inclination is fixed at i s = 90 ◦ while α is changed from 0.0 up to 0 .
4. The RM curvesfor various stellar inclinations are shown in the right panel along with the case for rigid rotation.Different colors indicate different coefficients ( α ) and stellar inclinations ( i s ) as shown in each panel.For each panel, RV residuals from the rigid rotation (the black line) are shown at the bottom.In order to evaluate the impact of differential rotations, we take XO-3 again as an example.Since the XO-3 system is reported to have a large spin-orbit misalignment and its host star has alarge rotational velocity, we can expect a large impact of the differential rotation on RM velocityanomalies. Figure 8 indicates the RM anomaly curves for XO-3 with and without differentialrotation. In the left panel, we change the differential rotation coefficient α , fixing the stellarinclination i s . We try the case of α = 0 .
0, 0.1, 0.2, 0.3, and 0.4. The right panel of Figure 8indicates the variation of the RM velocity anomaly for the different cases of i s , in which we set i s = 90 ◦ , 60 ◦ , and 30 ◦ while α is fixed at 0.2. In this figure, the stellar spin velocity at the equatoris fixed as v sin i s = R s Ω eq sin i s = 18 . − . At the bottom of each panel, we show the residualof each curve from the black line (rigid rotation).It is notable that the residuals are comparable to or rather larger than the RV precision for 27 –this system, which is reported to be ∼ − (Winn et al. 2009). The largest impact ( ∼
15 ms − ) on ∆ v occurs for the case of large values of α and small values of i s , as is expected. Figure 8also shows that the stellar inclination i s significantly affects the velocity anomaly. This result, inturn, suggests that we may be able to estimate the stellar inclination i s from RM analyses if wecan model the differential rotation for each type of stars and fix the coefficient α .
5. Discussion and Summary
We have developed a new analytic formula (Eq. [16]) to describe the velocity anomaly duringa planetary transit. We here summarize the previous findings describing the RM velocity anomalyalong with our new analytic formula:∆ v = − f − f v p , (OTS) (26)∆ v ≈ − (cid:26) β ⋆ β ⋆ + β p (cid:27) / f v p " − v p β ⋆ + β p + v p β ⋆ + β p ) , (Hirano et al. 2010) (27)∆ v ≈ − f π Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) ˜Θ( σ ) sin(2 πσv p ) σdσ Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) n ˜ M ( σ ) − f ˜Θ( σ ) cos(2 πσv p ) o σ dσ , (this work)where β ⋆ and β p in Equation (27) indicate the best-fit Gaussian dispersions for the stellar templatespectrum (including stellar rotation), and for the intrinsic line broadening in the absence of stellarrotation, respectively (see Hirano et al. (2010) for details). Major improvement in our new formulais that it precisely modeled various effects which can possibly affect stellar line profiles such asrotational broadening, macroturbulence, and the instrumental profile.Numerical simulations using the mock transit spectra and the RV analysis routine for Sub-aru/HDS have shown that our analytic formula agrees with the simulated results within ∼ .
5% oftypical RM amplitudes ( . .
25 m s − in the case of a typical transit of a hot Jupiter).We have also compared Equation (16) with the output of numerical calibrations of the RMeffect that have been given in the literature, based on mock data simulations for Keck/HIRES.We have plotted the RM velocity anomaly curves for four existing transiting systems: HAT-P-4,HAT-P-14, XO-3, and HAT-P-11. Equation (16) proved to be in good agreement with previouslyreported empirical relations, although some deviations were seen for the XO-3 system. We suggestthat this disagreement comes from the coupling between rotational and macroturbulent broadening,and perhaps also to imperfectly known stellar line parameters. One possible means by which thisissue could be clarified is to apply our new analytic formula to observed RVs and see if Equation(16) gives a better fit to the data. We do not attempt to do so in the present paper, but we intendto reanalyze many of the existing data-sets in due course.As Hirano et al. (2010) and others have pointed out, the relation between the position of the 28 –planet and the observed RM velocity anomaly usually does not significantly change the estimationfor the spin-orbit misalignment angle λ since a spin-orbit misalignment is most likely to be shown byan asymmetry in the velocity anomaly curve during a transit. However, when the transit impactparameter b is small, the asymmetry in the velocity anomaly curve becomes less clear and theamplitude of a series of ∆ v becomes more important.Analytic approaches to a phenomenon give us an insight into the behavior of the phenomenon,and sometimes reveals important aspects which numerical approaches are likely to overlook. Wehope that our new formula will be useful in data analyses of the RM effect and make a contributionto the progress in this field.We wish to acknowledge Dr. Paula R.T. Coelho for her kind instruction on the use of thesynthetic spectra. We are very grateful to John Asher Johnson for helpful discussions on thistopic. The data analysis was in part carried out on common use data analysis computer systemat the Astronomy Data Center, ADC, of the National Astronomical Observatory of Japan. T.H.is supported by Japan Society for Promotion of Science (JSPS) Fellowship for Research (DC1: 22-5935). Y.S. gratefully acknowledges support from the Global Collaborative Research Fund (GCRF)“A World-wide Investigation of Other Worlds” grant and the Global Scholars Program of PrincetonUniversity, and also from JSPS Core-to-Core Program “International Research Network for DarkEnergy”. J.N.W. acknowledges support from the NASA Origins program (NNX11AG85G). N.N.acknowledges a support by NINS Program for Cross-Disciplinary Study. We would like to expressspecial thanks the anonymous referee for the helpful comments and suggestions on this manuscript. A. Derivation of the Analytic Formula (16)
In this appendix, we derive Equation (16) based on Equations (1), (10), (11), and (12).In order to calculate convolutions and cross-correlations, dealing with them in the Fourierdomain, where they are expressed as products, can significantly facilitate the computation. First,we Fourier-transform Equations (1) and (10) as˜ F star ( σ ) = − ˜ S ( σ ) ˜ M ( σ ) , (A1)˜ F transit ( σ ) = − ˜ S ( σ ) ˜ M ( σ ) + f ˜ S ( σ ) ˜Θ( σ ) e πiσv p , (A2)where the Fourier transform of an arbitrary function F ( v ) is defined as˜ F ( σ ) ≡ Z ∞−∞ F ( v ) e − πiσv dv, (A3)so that ˜ M ( σ ) and ˜Θ( σ ) are explicitly given by Equations (17) and (18), respectively (see Appendix 29 –B). Thus, the Fourier transform of the cross-correlation function C ( x ) becomes˜ C ( σ ) = ˜ F star ( σ ) · ˜ F ∗ transit ( σ )= n ˜ S ( σ ) o ˜ M ( σ ) h ˜ M ( σ ) − f ˜Θ( σ ) e − πiσv p i . (A4)The Fourier transforms of the Gaussian and Lorentzian functions respectively are written as G ( v ) = 1 β √ π e − v /β = ⇒ ˜ G ( σ ) = e − π β σ , (A5) L ( v ) = 1 π γ v + γ = ⇒ ˜ L ( σ ) = e − πγ | σ | . (A6)When we assume the intrinsic line profile as S ( v ) = V ( v ; β, γ ) = G ( v : β ) ∗ L ( v ; γ ) (the Voigtfunction), the Fourier transformation of S ( v ) becomes˜ S ( σ ) = e − π β σ − πγ | σ | . (A7)In this case, the inverse Fourier transformation for ˜ C ( σ ) is expressed as C ( x ) = Z ∞−∞ ˜ C ( σ ) e πiσx dσ = Z ∞−∞ exp( − π β σ − πγ | σ | + 2 πiσx ) ˜ M ( σ ) (cid:16) ˜ M ( σ ) − f ˜Θ( σ ) e − πiσv p (cid:17) dσ = 2 Z ∞ exp( − π β σ − πγσ ) n ˜ M ( σ ) o cos(2 πσx ) dσ − f Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) ˜Θ( σ ) cos { πσ ( v p − x ) } dσ. (A8)The derivative of Equation (A8) with respect to x is dC ( x ) dx = − π Z ∞ exp( − π β σ − πγσ ) n ˜ M ( σ ) o σ sin(2 πσx ) dσ − πf Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) ˜Θ( σ ) σ sin { πσ ( v p − x ) } dσ. (A9)The maximum of the cross-correlation function is at x = ∆ v , where Z ∞ exp( − π β σ − πγσ ) n ˜ M ( σ ) o σ sin(2 πσ ∆ v ) dσ = − f Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) ˜Θ( σ ) σ sin { πσ ( v p − ∆ v ) } dσ. (A10)Because of the exponential factors, the integrals in Equation (A10) can be cut off at σ ∼ / β .Moreover, since the Gaussian width parameter β is of the order of several km s − and the velocityanomaly is typically less than one hundred m s − , we obtain2 πσ ∆ v . . , (A11) 30 –and we can safely approximate as sin(2 πσ ∆ v ) ≈ πσ ∆ v . Thus, using the identity for trigonometricfunctions: sin { πσ ( v p − ∆ v ) } = sin(2 πσv p ) cos(2 πσ ∆ v ) − cos(2 πσv p ) sin(2 πσ ∆ v ) ≈ sin(2 πσv p ) − πσ ∆ v cos(2 πσv p ) , (A12)we obtain following analytic formula for ∆ v :∆ v ≈ − f π Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) ˜Θ( σ ) sin(2 πσv p ) σdσ Z ∞ exp( − π β σ − πγσ ) ˜ M ( σ ) n ˜ M ( σ ) − f ˜Θ( σ ) cos(2 πσv p ) o σ dσ . (A13)Note that contrary to previously reported empirical relations, Equation (A13) implies that thevelocity anomaly due to the RM effect is not proportional to the loss of flux f , but the secondorder of f slightly affects the results. B. The Fourier Transformation of the Rotational and Macroturbulence Kernel
One of the great successes in dealing with the RM effect in the Fourier domain is that therotational and macroturbulence broadening kernel M ( v ) becomes significantly simpler in the Fourierdomain. Here, we Fourier-transform M ( v ) and derive Equation (17).Since the velocity component v appears only in the macroturbulence kernel Θ( v ), we can writethe Fourier transformation of M ( v ) as˜ M ( σ ) ≡ Z ∞−∞ M ( v ) e − πiσv dv ∝ Z ∞−∞ Θ( v − x Ω sin i s ) e − πiσv dv. (B1)The last expression above is expressed by a linear combination of the following form of an integral: Z ∞−∞ √ πκ exp ( − (cid:18) v − x Ω sin i s κ (cid:19) − πiσv ) dv, (B2)where κ is a constant. The exponent in the above integrand reduces to − (cid:18) v − x Ω sin i s κ (cid:19) − πiσv = − ( v − x Ω sin i s + πiκ σ ) κ + πσi ( πκ σi − x Ω sin i s ) . (B3)Thus, using the Gaussian integral: Z ∞−∞ √ πκ e − ( v/κ ) dv = 1 , (B4) 31 –Equation (B2) becomes Z ∞−∞ √ πκ exp ( − (cid:18) v − x Ω sin i s κ (cid:19) − πiσv ) dv = exp { πiσ ( πiκ σ − x Ω sin i s ) } . (B5)Replacing κ with ζ cos θ or ζ sin θ , the second line in Equation (B1) is expressed as Z ∞−∞ Θ( v − x Ω sin i s ) e − πiσv dv = e − πiσx Ω sin i s { e − ( πζ cos θ ) σ + e − ( πζ sin θ ) σ } . (B6)Therefore, we obtain the following expression for ˜ M ( σ ):˜ M ( σ ) = Z Z entire disk − u (1 − cos θ ) − u (1 − cos θ ) π (1 − u / − u / × e − πiσx Ω sin i s { e − ( πζ cos θ ) σ + e − ( πζ sin θ ) σ } dxdyR s . (B7)Since the imaginary part of the integrand in Equation (B7) is an odd function with regard to x ,it vanishes after integrated over x . We then change variables as x = r cos φ and y = r sin φ , whichleads to ˜ M ( σ ) = Z R s Z π − u (1 − cos θ ) − u (1 − cos θ ) π (1 − u / − u / × cos(2 πσr Ω sin i s cos φ ) e − ( πζ cos θ ) σ + e − ( πζ sin θ ) σ rdrdφR s , (B8)where cos θ = s − r R s , sin θ = rR s . (B9)Using the following formula Z π cos(2 πσr Ω sin i s cos φ ) dφ = 2 πJ (2 πσr Ω sin i s ) , (B10)we obtain ˜ M ( σ ) = Z R s − u (1 − cos θ ) − u (1 − cos θ ) − u / − u / × { e − ( πζ cos θ ) σ + e − ( πζ sin θ ) σ } J (2 πσr Ω sin i s ) rdrR s , (B11)or equivalently, ˜ M ( σ ) = Z − u (1 − √ − t ) − u (1 − √ − t ) − u / − u / × { e − π ζ σ (1 − t ) + e − π ζ σ t } J (2 πσv sin i s t ) tdt. (B12) 32 – C. Estimation of the Effective Line Profile
Here, we describe how to estimate the “effective” line profile from a spectrum. First, wemodel the auto-correlation function of a spectral assuming a simple analytic form of it, and thencompute the actual auto-correlation function for the synthetic spectra that we use in the mock datasimulation.We assume that the synthetic spectrum (Coelho et al. 2005) is expressed by a summation ofVoigt functions (defined by Eq. [13]) as T ( λ ) = X i k i V ( λ − d i ; β i , γ i ) , (C1)where k i , d i are the equivalent width and the central wavelength of the spectral line i . The Gaussianand Lorentzian dispersions of the line i are also noted by β i and γ i , respectively. For simplicity, wesubtract the continuum (normalized to unity) of the spectrum and flip the absorption lines so thatthe flux takes positive values. The auto-correlation function A ( x ) of the spectrum is defined as A ( x ) ≡ Z ∞−∞ T ( λ ) T ( λ − x ) dλ = Z ∞−∞ X i X j k i k j V ( λ − d i ; β i , γ i ) V ( λ − d j − x ; β j , γ j ) dλ. (C2)Since the convolution between two Voigt functions yields another Voigt function with a largerdispersion, the integral above reduces to A ( x ) = X i X j k i k j V ( x − d i + d j ; q β i + β j , γ i + γ j )= X i k i V ( x ; √ β i , γ i ) + X i X i>j k i k j V ( x − | d i − d j | ; q β i + β j , γ i + γ j )+ X i X i>j k i k j V ( x + | d i − d j | ; q β i + β j , γ i + γ j ) . (C3)The first term in the last expression indicates the sum of all the absorption lines broadened byfactors of √ β and ¯ γ as X i k i V ( x ; √ β, γ ) ≡ X i k i V ( x ; √ β i , γ i ) . (C4) 33 –Although a summation of Voigt functions is not necessarily expressed by a Voigt function, thistreatment above is justified as long as all the values of β i and γ i do not differ significantly from lineto line.The second and third terms in Equation (C3) mean the correlations between the lines separatedby | d i − d j | . It is plausible to assume that this separation | d i − d j | is distributed randomly in aspectrum. In this case, the second and third terms become the summations of a bunch of Voigtfunctions centered at randomly placed x . Assuming the number of spectral lines is large enough( i & x .Therefore, we can approximate Equation (C3) around x = 0 as A ( x ) ≈ ¯ kV ( x ; √ β, γ ) + C, ¯ k ≡ X i k i , (C5)where C is the constant originated from the second and third terms in Equation (C3).Next, we try to compute the auto-correlation function A ( x ) for the actual synthetic spectrawe use to generate the mock transit spectra. We take the G0-type spectrum described in Section3.1.1, subtract the continuum, and flip the lines. Then, we compute the auto-correlation function A ( x ) defined as Equation (C2). Note that we cut the spectrum so that it ranges from ∼ ∼ A ( x ) x [Angstrom] Fig. 9.— The auto-correlation function A ( x ) of the synthetic spectrum for a G0-type star, definedby Equation (C2). The wavelength range used in the calculation is 5000 ˚A . λ . A ( x ) of the synthetic spectrum for a G0- In reality, since we deal with a limited wavelength range of the spectrum, the auto-correlation function A ( x )slowly decreases as a function of | x | .
34 –type star used in the mock data simulation around x = 0. As is expected from the above theoreticalprediction, it has a strong Voigt-shaped peak at x = 0 with an offset due to the correlations withneighboring lines. We fit this auto-correlation function assuming a form of Equation (C5). Sincethe Gaussian and Lorentzian dispersion β and γ are strongly correlated, we fix ¯ β as ¯ β = β = 1 . − , which is calculated by Equation (19), and let only ¯ k , ¯ γ , and C be free. As a result, wefind the best-fit value of ¯ γ as ¯ γ ∼ .
94 km s − . In this paper, we use these values in comparing theanalytic formula with the simulated results for a G0-type star. REFERENCES
Albrecht, S., Reffert, S., Snellen, I., Quirrenbach, A., & Mitchell, D. S. 2007, A&A, 474, 565Butler, R. P., Marcy, G. W., Williams, E., McCarthy, C., Dosanjh, P., & Vogt, S. S. 1996, PASP,108, 500Chatterjee, A., Sarkar, A., Barat, P., Mukherjee, P., & Gayathri, N. 2008, ApJ, 686, 580Claret A. 2004, A&A, 428, 1001Coelho, P., Barbuy, B., Mel´endez, J., Schiavon, R. P., & Castilho, B. V. 2005, A&A, 443, 735Collier Cameron, A., Guenther, E., Smalley, B., McDonald, I., Hebb, L., Andersen, J., Augusteijn,Th., Barros, S. C. C., Brown, D. J. A., Cochran, W. D., Endl, M., Fossey, S. J., Hartmann,M., Maxted, P. F. L., Pollacco, D., Skillen, I., Telting, J., Waldmann, I. P., & West, R. G.2010, MNRAS, 407, 507Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298Fabrycky, D. C., & Winn, J. N. 2009, ApJ, 696, 1230Gaudi, B. S., & Winn, J. N. 2007, ApJ, 655, 550Gim´enez, A. 2006, ApJ, 650, 408Gray, D. F. 2005, The Observation and Analysis of Stellar Photospheres (3rd ed.; CambridgeUniversity Press.)Hirano, T., Suto, Y., Taruya, A., Narita, N., Sato, B., Johnson, J. A., & Winn, J. N. 2010, ApJ,709, 458Hirano, T., Narita, N., Shporer, A., Sato, B., Aoki, W., & Tamura, M. 2011, PASJ, 63, S531Hosokawa, Y. 1953, PASJ, 5, 88Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606 35 –Lubow, S. H., & Ida, S. 2010, arXiv:1004.4137McLaughlin, D. B. 1924, ApJ, 60, 22Morton, T. D., & Johnson, J. A. 2011, ApJ, 729, 138Moutou, C., Diaz, R. F., Udry, S., Hebrard, G., Bouchy, F., Santerne, A., Ehrenreich, D., Arnold,L., Boisse, I., Bonfils, X., Delfosse, X., Eggenberger, A., Forveille, T., Lagrange, A., Lovis,C., Martinez, P., Pepe, F., Perrier, C., Queloz, D., Santos, N. C., Segransan, D., Toublanc,D., Troncin, J., Vanhuysse, M., & Vidal-Madjar, A. 2011, A&A, in pressNagasawa, M., Ida, S., & Bessho, T. 2008, ApJ, 678, 498Narita, N., Enya, K., Sato, B., Ohta, Y., Winn, J. N., Suto, Y., Taruya, A., Turner, E. L., Aoki,W., Yoshii, M., Yamada, T., & Tamura, Y. 2007, PASJ, 59, 763Narita, N., Hirano T., Sato B., Winn, J. N., Suto, Y., Turner, E. L., Aoki, W., Tamura, M., &Yamada, T. 2009a, PASJ, 61, 991Narita, N., Sato, B., Hirano, T., & Tamura, M. 2009b, PASJ, 61, L35Narita, N., Sato, B., Hirano, T., Winn, J. N., Aoki, W., & Tamura, M. 2010, PASJ, 62, 653Narita, N., Hirano, T., Sanchis-Ojeda, R., Winn, J. N., Holman, M. J., Sato, B., Aoki, W., &Tamura, M. 2010b, PASJ, 62, L61Ohta, Y., Taruya, A., & Suto Y. 2005, ApJ, 622, 1118Queloz, D., Eggenberger, A., Mayor, M., Perrier, C., Beuzit, J. L., Naef, D., Sivan, J. P., & Udry,S. 2000, A&A, 359, L13Reiners, A. 2003, A&A, 408, 707Rossiter, R. A. 1924, ApJ, 60, 15Sato, B., Kambe, E., Takeda, Y., Izumiura, H., & Ando, H. 2002, PASJ, 54, 873Shporer, A., & Brown, T. 2011, 733, 30Valenti, J. A., & Fischer, D. A. 2005, ApJ Supplement Series, 159, 141Triaud, A. H. M. J., et al. 2009, A&A, 506, 377Winn, J. N., Noyes, R. W., Holman, M. J., Charbonneau, D., Ohta, Y., Taruya, A., Suto, Y.,Narita, N., Turner, E. L., Johnson, J. A., Marcy, G. W., Butler, R. P., & Vogt S. S. 2005,ApJ, 631, 1215Winn, J. N., Johnson, J. A., Fabrycky, D., Howard, A. W., Marcy, G. W., Narita, N., Crossfield,I. J., Suto, Y., Turner, E. L., Esquerdo, G., & Holman, M. J. 2009, ApJ, 700, 302 36 –Winn, J. N. 2010, arXiv:1001.2010Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010a, ApJL, 718, L145Winn, J. N., Johnson, J. A., Howard, A. W., Marcy, G. W., Isaacson, H., Shporer, A., Bakos, G.´A., Hartman, J. D., & Albrecht, S. 2010b, ApJL, 723, L223Winn, J. N., Howard, A. W., Johnson, J. A., Marcy, G. W., Isaacson, H., Shporer, A., Bakos, G.´A., Hartman, J. D., Holman, M. J., Albrecht, S., Crepp, J. R., & Morton, T. D. 2011, ApJ,141, 63Wu, Y., Murray, N. W., & Ramsahai, J. M. 2007, ApJ, 670, 820