TRANSIT LIGHTCURVES OF EXTRASOLAR PLANETS ORBITING RAPIDLY ROTATING STARS

In order to generate synthetic transit lightcurves with fast-rotating stars, I modified the algorithm originally developed for Barnes & Fortney (2003) and extended in Barnes & Fortney (2004) and Barnes (2007). The algorithm numerically integrates the total flux coming from the uneclipsed star, F₀, in polar coordinates centered on the projected center of the star in the plane of the sky such that

$$\begin{equation} F_0 = \int _{0}^{R_\mathrm{eq}} \int _{0}^{2\pi } I(r, \theta) d\theta dr, \end{equation} \tag{ 1 }$$

where R_eq is the radius of the star at its equator (see Figure 1 for a schematic of some of the geometric variables), r and θ are measured from the stellar center and counterclockwise from the x-axis respectively, and I(r, θ) is the star's intensity at point (r, θ). It then evaluates the apparent stellar flux at time t, F(t), relative to the out-of-transit flux F₀, by subtracting the amount of stellar flux blocked by the planet from F₀:

$$\begin{equation} F_{\rm blocked}(t) = \int _0^{R_\mathrm{eq}} \int _{0}^{2\pi } \Gamma (r,\theta,t) I(r,\theta) d\theta dr \mathrm{} \end{equation} \tag{ 2 }$$

and

$$\begin{equation} F(t) = \frac{F_0 - F_{\rm blocked}}{F_0} \mathrm{,} \end{equation} \tag{ 3 }$$

where Γ(r, θ, t) = 1 if the planet is blocking starlight at position r, θ, and time t, and Γ(r, θ, t) = 0 if not.

Zoom In Zoom Out Reset image size

The difference in the case of fast-rotating stars is that the rotation induces these stars' equators to bulge outward. Hence the polar integral in Equations (1) and (2) no longer properly accounts for the symmetry of the problem. I quantify this effect here as the star's oblateness, f, defined to be $f\equiv \frac{R_\mathrm{eq}-R_\mathrm{pole}}{R_\mathrm{eq}}$ where R_pole is the star's radius along the rotational pole. I assume throughout that the star's resulting shape can be considered as a MacLaurin spheroid. The value relevant for the integrals in Equations (1) and (2), though, is the effective oblateness f_eff, which I define to be the apparent oblateness of the star when projected into the plane of the sky. The effective oblateness is related to the actual oblateness by the stellar obliquity φ, where φ = 0 if the stellar rotation axis resides in the plane of the sky:

$$\begin{equation} f_\mathrm{eff} = 1 - \sqrt{(1-f)^2\cos ^2\varphi + \sin ^2\varphi }. \end{equation} \tag{ 4 }$$

The relationship is not a simple cosine owing to three-dimensional geometry. I will show how this expression can be derived a few paragraphs below.

With f_eff in hand, it becomes straightforward to incorporate the stellar asphericity. In order to avoid complex and computationally intensive elliptical integrals, I execute a substitution for r and θ in Equations (1) and (2). I instead choose to integrate over r' and θ', where r' and θ' are chosen so as to "pop" the star into a spherical shape in r'–θ' space. To do this, I first convert the true projected r and θ measured from the star's center into x = rcos θ and y = rsin θ. I assume that the projected stellar rotation axis is parallel to the y-axis for simplicity—the true orientation will not be known, in general, but does not matter since it does not affect the measured stellar flux. I then let x' = x and

$$\begin{equation} y^{\prime }=\frac{y}{(1-f_\mathrm{eff})}, \end{equation} \tag{ 5 }$$

and then set $r^{\prime }=\sqrt{x^{\prime 2}+y^{\prime 2}}$ and θ' = atan2(y', x'), where atan2 is the C computer language arctangent function that returns a true 4-quadrant-capable angle from x and y values. The substitution in Equation (5) works equally well if you were to choose to integrate the stellar flux in Cartesian xy-coordinates instead of the polar integral that I use.

Equations (1) and (2) now become

$$\begin{equation} F_0 = (1-f_\mathrm{eff}) \int _{0}^{R_\mathrm{eq}} \int _{0}^{2\pi } I(r^{\prime }, \theta ^{\prime }) d\theta ^{\prime } dr^{\prime } \end{equation} \tag{ 6 }$$

and

$$\begin{equation} F_{\rm blocked}(t) = (1-f_\mathrm{eff}) \int _0^{R_\mathrm{eq}} \int _{0}^{2\pi } \Gamma (r^{\prime },\theta ^{\prime },t) I(r^{\prime },\theta ^{\prime }) d\theta ^{\prime } dr^{\prime } \mathrm{.} \end{equation} \tag{ 7 }$$

The (1 − f_eff) factor is introduced by the coordinate transformation: "popping" the oblate star out into a circle overestimates its projected area, and hence the emitted flux, by (1 − f_eff)⁻¹. In the end the factor is irrelevant. When the results are plugged into Equation (3), it drops out. Hence the lightcurve generation algorithm as implemented does not use the factor (1 − f_eff) explicitly at all.

All that is left then is to determine I(r', θ') and then to integrate it. Obtaining I(r, θ) is straightforward but nontrivial. I first break out the stellar limb darkening from the normal emission and assume blackbody radiation

$$\begin{equation} I(r^{\prime },\theta ^{\prime }) = B_\lambda (T(r^{\prime },\theta ^{\prime })) L(r^{\prime }, \theta ^{\prime }) \end{equation} \tag{ 8 }$$

where B_λ(T) is the blackbody function (or your desired stellar emission as a function of temperature at a given wavelength), T(r', θ') is the temperature at a given point on the stellar disk, and L(r', θ') is the stellar limb darkening at that point. It would be possible to incorporate a more realistic stellar emission spectrum rather than to assume it to be a blackbody, but since rapidly rotating stars are mostly of early spectral types, the differences are not significant at the level of the present investigation.

A modified version of the von Zeipel theorem (von Zeipel 1924) determines the stellar temperature T at every point. The temperature on the surface of a rapidly rotating star is given by Maeder (2009):

$$\begin{equation} T = T_\mathrm{pole} \frac{g^\beta }{g_\mathrm{pole}^\beta }, \end{equation} \tag{ 9 }$$

where g is the magnitude of the local effective surface gravity, and T_pole and g_pole are the pole's temperature and surface gravity, respectively. The value β is known as the gravity-darkening parameter. Its nominal value is 0.25 for a purely radiative star, as derived by von Zeipel; for our numerical calculations I use an empirically measured β, as detailed below. The local surface gravity vector $\vec{g}$ has two terms: one Newtonian and one centrifugal:

$$\begin{equation} \vec{g} = -\frac{GM_*}{R^2}\hat{r} + \Omega ^2 R_\perp \hat{r_\perp }, \end{equation} \tag{ 10 }$$

where G is the universal gravitational constant, M_* is the stellar mass, Ω is the stellar rotation rate in radians per second, R is the distance from the star center to the point in question, and R_⊥ is the distance from the star's rotation axis to the point in question. The two hatted symbols, $\hat{r}$ and $\hat{r_\perp }$ are unit vectors pointing to the point in question from the stellar center and the stellar rotation axis, respectively.

So, then, to know I you need to know g, and to know g you need to know the three-dimensional vector position of each point that you see on the star. In this case, you already know r' and θ', from which you can trivially derive x and y, the location of two-dimensional projection in the plane of the sky of the point of interest with respect to the center of the star. What is left then is to determine z. This is nontrivial.

The geometrical constraint on z is that its value must conform to the surface of an oblate spheroid given x, y, the stellar radius R_eq, and the oblateness f. This is easy enough when the stellar obliquity is zero, and thus when the y-axis is parallel to the stellar rotation axis. So I define a new set of coordinate axes with the same origin as the x–y–z system, at the center of the star. However, this new system is rotated in the y–z plane (i.e., around the x-axis) by an angle φ, the star's obliquity to the plane of the sky. Call this the x₀, y₀, z₀ system, where

$$\begin{eqnarray} x_0 & \equiv & x,\nonumber \\ y_0 & \equiv & y\cos (\varphi) + z\sin (\varphi), \mathrm{and}\nonumber \\ z_0 & \equiv & -\! y\sin (\varphi) + z\cos (\varphi). \end{eqnarray} \tag{ 11 }$$

In this new obliquity rotated coordinate system, the surface of the star's photosphere follows

$$\begin{equation} x_0^2+\frac{y_0^2}{(1-f)^2}+z_0^2=R_\mathrm{eq}^2. \end{equation} \tag{ 12 }$$

Plugging the definitions of x₀, y₀, and z₀ from Equations (11) into Equation (12), I solve for z in terms of known x and y. The solution to the resulting quadratic is

$$\begin{equation} z = \frac{-2y(1-(1-f)^2)\sin {\varphi }\cos {\varphi }+\sqrt{d}}{2((1-f)^2\cos ^2\varphi +\sin ^2\varphi)}, \end{equation} \tag{ 13 }$$

where the determinant d is

$$\begin{eqnarray} \hspace{-25pt}d &\equiv & 4y^2(1-(1-f^2))^2\sin ^2\varphi \cos ^2\varphi \nonumber\\ \hspace{-25pt}&&-4((\cos ^2\varphi (1-f)^2+\sin ^2\varphi) \nonumber \\ \hspace{-25pt}&&\times\big(\big(y^2\sin ^2\varphi -R_\mathrm{eq}^2+x^2\big)(1-f^2)+y^2\cos ^2\varphi)\big). \end{eqnarray} \tag{ 14 }$$

I choose the positive root of the determinant as the negative root represents the invisible second interception of the line of sight with the photosphere that occurs on the backside of the star as seen from the Earth. Equation (4) can be derived by setting the determinant d equal to zero to establish the outer edge of the star's disk as seen from the Earth, and then solving for the proper f_eff to reproduce that disk.

Now I have all of the necessary parameters to compute the flux coming from each point on the star. To do so, at each x, y point use Equation (13) to get z, and then plug z into Equations (11) to get x₀, y₀, and z₀. Do the vector addition to get $\vec{g}$ from Equation (10) in x₀, y₀, z₀ space where

$$\begin{eqnarray*} \vec{R} &\equiv &\frac{x_0}{R}\hat{i_0}+\frac{y_0}{R}\hat{j_0}+\frac{z_0}{R}\hat{k_0},\\ \vec{R_\perp } &\equiv & \frac{x_0}{R_\perp }\hat{i_0}+\frac{z_0}{R_\perp }\hat{k_0},\\ R & \equiv & \sqrt{x_0^2+y_0^2+z_0^2}, \\ R_\perp & \equiv & \sqrt{x_0^2+z_0^2}. \end{eqnarray*}$$

Then plug $g\equiv |\vec{g}|$ into von Zeipel's equation (Equation (9)) to get T, and then derive a flux from T using a blackbody curve or your choice of a more sophisticated emitted flux.

In order to generate appropriate and representative lightcurves that can be used for comparison to transits yet undiscovered, I calculate all transit lightcurves as if the parent star were Altair (α Aquilae). The true host stars for Kepler-detected transiting planets will show varying stellar masses, polar temperatures, radii, and rotation rates. I elect to use Altair because I think that its spectral type (A7V) is broadly representative of the majority of the expected fast-rotating stars in the Kepler sample, and because its parameters are well characterized by interferometric imaging (Monnier et al. 2007). Not all sets of parameters produce physically plausible stars; I avoid non-physical combinations by only using this one set of known stellar values.

The specific stellar parameters that I use are M_* = 1.8 M_☉ (Peterson et al. 2006a), R_eq = 2.029 R_☉, T_pole = 8450 K, β = 0.190, f = 0.1947, and a stellar rotation period of 8.64 hr (all as directly measured for Altair by Monnier et al. 2007). For the planet, I assume R_p = R_Jup and an orbit semimajor axis of 0.05 AU (corresponding to a period of 3.04 days) for familiarity with the lightcurves of known transiting hot Jupiters. Transit lightcurve shapes are invariant with orbit period; only the timescale changes. Hence the curves that I show here can be converted for different-period planets by stretching the x-axis.

In fitting I adjust for the c₁ and c₂ limb darkening coefficients outlined in Brown et al. (2001). I generate the synthetic lightcurves using c₁ = 0.640 and c₂ = 0.0. I also assume monochromatic observation at 0.51 μm wavelength except where otherwise noted.

Just because the transit lightcurve of a planet around a fast-rotating star is symmetric does not mean that the planet is spin–orbit aligned. Any planet whose transit chord is perpendicular to the projected stellar rotation pole, i.e. for which the angle between the transit chord and the projected pole, α, is 00, will show a symmetric lightcurve. I plot the lightcurve for a few such hypothetical planets in Figures 3 and 4. Figure 3 shows planets around a star with obliquity φ = 30°. Figure 4 shows transits for a star viewed pole-on with φ = 90°.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although these lightcurves are symmetrical about the mid-transit point, unlike the spin-aligned case they are not well-fit using a model that assumes a spherical star. In the φ = 30° case the stellar obliquity removes the symmetry in the stellar disk, thereby making transits toward the north and south poles distinct from one another. I arbitrarily choose to define those transits toward the north of the center of the stellar disk to have a negative impact parameter so as to differentiate the two possible cases.

In the positive impact parameter cases in Figure 3, the planet crosses the cooler equatorial regions of the star. Hence those transits are relatively shallow. The geometry of the orthographic projection of the star as seen from Earth means that for the high impact parameter, b = 0.6 R_pole and b = 0.9 R_pole cases, the planet covers more southerly stellar latitudes at ingress and egress than it does at mid-transit. Because those more southerly latitudes are hotter, the net result is a flatter lightcurve bottom than occurs for spherical stars with the same limb darkening parameter (c₁ = 0.640). The flatness appears as a lower limb darkening parameter c₁ in Table 1. For the b = 0.6 R_pole and b = 0.9 R_pole cases the model overestimates the star's radius; the best-fit radius is greater than even the star's equatorial radius.

For b = 0.0 R_pole and b = 0.3 R_pole in the φ = 30° case, the planet transits closer to the star's north pole, and hence the transit depth becomes progressively greater. The ingress and egress are again at more southerly latitudes than mid-transit. But now mid-transit is north of the equator, so the ingress and egress are cooler than mid-transit, enhancing the curvature of the lightcurve bottom. I was unable to find satisfactory fits to either of these synthetic lightcurves using just the first Brown et al. (2001) limb darkening coefficient c₁ alone, even though it was the only one used to generate the synthetic lightcurve. However, by allowing the fitting algorithm to optimize c₂ as well I managed to find reasonable fits for b = 0.0 R_pole and b = 0.3 R_pole.

When fitting b = −0.3 R_pole, b = −0.6 R_pole, and b = −0.9 R_pole not even two-limb-darkening fits were acceptable. The curvature on these transits is so extreme that even though they are symmetric, they are discernable from reasonable spherical-star models because of the severity of their curvature. In fact, the curvature of these planets' transit bottoms are so severe that the lightcurves look "V"-shaped instead of "U"-shaped. The transit bottom is difficult to discern. The resulting rounded shape might be mistaken for an eclipsing binary star.

Fits for transits of pole-on oriented fast-rotating stars are similar. As shown in Figure 4, the center of the stellar disk is particularly bright relative to the limb owing to the center being the hot pole and the limb being the cool equator. Thus here again the curvature of the transit lightcurve between second and third contacts is large, and cannot be satisfactorily fit by a spherical-star model. In both these cases the difficulty in identifying these type of planetary transits in the Kepler data will probably be recognizing that they are not eclipsing binaries, not in thinking that they are planets around spherical stars.

An observer would also measure a symmetrical transit lightcurve if a planet orbits a zero obliquity (φ = 00) star such that its orbit plane contains the stellar rotation pole, i.e., with α = 90°. I show synthetic lightcurves for such a situation in Figure 5. These are strange transits. By moving parallel to the star's orbit pole, the planets first encounter the limb-darkened edge of the star. Shortly after second contact, they cover the brightest part of the stellar photosphere underneath the transit chord. The lightcurve is deepest there. Over the equator at mid-transit the transit depth is shallow, and then the process repeats itself backwards on egress.

Zoom In Zoom Out Reset image size

The resulting lightcurve has a "double-horned" structure that bears a resemblance to what a transit across a limb-brightened star would look like. However, spherical-star model fits cannot reproduce the specific structure. In particular, a limb-brightened star would have sharp points at second and third contacts where the derivative is discontinuous and changes sign. The fast-rotating star φ = 0°, α = 90° transits in Figure 5 instead have rounded horns at second and third contact owing to the combination of limb- and gravity-darkening. The resulting rounded double-horn shape is both characteristic and diagnostic of this kind of transit.

If the stellar orientation and transit chord have more exotic geometries, then highly unusual asymmetric lightcurves result. A star with obliquity φ = 30° and transit chord azimuth α = 90° (motion parallel to the projected stellar rotation axis) presents an easy such case to visualize (Figure 6). Because the stellar disk presents only a left-right symmetry, these transits pass over very different photospheric temperatures in the first half of the transit relative to the second half.

Zoom In Zoom Out Reset image size

In the φ = 30°, α = 90°, impact parameter b = 0.0 R_pole case, the planet first crosses a darkened stellar limb before occulting the star's hot north pole. The resulting lightcurve is deepest at this point. After mid-transit the planet occults the cool equator for a shallower transit depth, but then the depth increases again near third contact as the planet covers higher southern latitudes near the southern stellar limb. The same process occurs in a more abbreviated fashion as the impact parameter b increases.

If the transit chord azimuth is oblique, then the lightcurves can become quite complex as shown in Figure 7. The central b = 0.0 R_pole transit resembles those from Figure 6. The more northerly transits (those defined to have negative impact parameters) also show strong lightcurve asymmetries, with deep first halves and shallowing second halves. Some of the lightcurves turn over again near third contact, but the ones with more negative impact parameters do not.

Zoom In Zoom Out Reset image size

For positive, more southerly impact parameters, the depth asymmetry decreases. Further from the hot pole, the temperatures under the transit chord are more uniform. The b = 0.9 R_pole transit is particularly interesting. It shows a nearly uniform depth in time, but while the ingress is long and gently curving, the egress is abrupt. A higher photometric precision would be required to definitively identify such a transit lightcurve as being one from a fast-rotating star than would be necessary for some of the more spectacularly asymmetric lightcurves.

With oblique azimuths also come uncentered lightcurves. Because of the stellar oblateness, the center of these transits does not occur at the time of the planet's inferior conjunction, that is, when the planet is closest to the Earth. The total discrepancy can be up to a few tenths of the total transit duration, depending on the stellar oblateness and the transit geometry.

The oblique azimuth transits of fast-rotating stars also introduce an asymmetry in the duration of transit ingress and egress. This effect comes about because the angle between the stellar limb and the transit chord is different in ingress and egress, thereby causing the transiting planet to take different amounts of time to cross the limb in each case.

Lightcurve asymmetries can occur in zero-obliquity (φ = 0°) stars as well when the transit chord azimuth is oblique. I show such a case in Figure 8, with φ = 0° and α = 60°. The central b = 0.0 R_pole transit is symmetric, similar to those from Figure 5. Non-central transits are asymmetric and resemble those from Figure 6 with one-half of the transit being deeper than the other. In this φ = 0° oblique case, the transits with negative impact parameter are the time inverse of those with positive impact parameter.

Zoom In Zoom Out Reset image size