DISCOVERY AND MASS MEASUREMENTS OF A COLD, 10 EARTH MASS PLANET AND ITS HOST STAR

For a period of just under three days on 2009 October 10–12, we were able to obtain observations using the High Resolution Instrument (HRI) visible imager on the EPOXI (Christiansen et al. 2011) spacecraft when it was located ∼0.1 AU from Earth. Observations with the EPOXI HRI were requested in an attempt to constrain the microlensing parallax effect and obtain precise mass measurements of the MOA-2009-BLG-266Lb planet and its host star. We could obtain these observations because our target field was able to provide a better test of newly installed pointing control software than a less dense stellar field.

The EPOXI data consist of 4127 50 s exposures with the "clear-6" filter. To minimize data transfer requirements, the data were taken as 128 × 128 and 256 × 256 pixel sub-frames. The first 3375 exposures used 128 × 128 sub-frames, and the last 752 images were 256 × 256 sub-frames. However, the pointing stability was such that the target occasionally drifted out of the 128 × 128 sub-frames, and it was only possible to do photometry on 2900 of these 3375 128 × 128 pixel images. An example of one of these 128 × 128 pixel images is shown on the right side of Figure 2.

The instrumental point-spread function (PSF) of the HRI on the Deep Impact probe is strongly dependent on the color of the target star. In addition, the instrument is permanently defocused, yielding a toroidal-shaped PSF. This is clearly not optimal for crowded field photometry, which typically requires a spatially varying empirical PSF model for deblending. Since we performed photometry using the Daophot/Allstar/Allframe software suite (Stetson 1994), we first required a model of the instrumental PSF that is usable by Daophot.

Instrumental PSFs have been generated by Barry (2010) for the HRI using the Drizzle algorithm (Fruchter & Hook 2002). In this process several hundred images of a standard star were added together to create a ten-time oversampled PSF model appropriate for that object. To approximate the PSF of MOA-2009-BLG-266, with V − I = 1.82, we co-added the instrumental PSFs of GJ436 with V − I = 2.44; and XO-2 with V − I = 0.75 with weightings of 0.795 and 0.205, respectively. This composite PSF was added to an otherwise empty image in a 3 × 3 grid, with each realization downsampled to standard resolution using a center pixel shifted by ±5 pixels in x and/or y. Daophot was then run on this image, using all nine images to build its own internal, double-resolution PSF model using an empirical function plus "lookup table."

Approximately 1% of pixels in our 50 s observations contain signatures of cosmic rays. We filtered these pixels using an algorithm that identifies features sharper than expected from the PSF, through pairwise comparison of neighboring pixels. These pixels were masked and objects underneath these pixels ignored when generating light curves. We used Daophot to detect stars in each image, and then Allstar to perform joint PSF photometry on all stars in a given image. We generated a master starlist by matching the results of the Allstar analysis using Daomatch and Daomaster. This starlist was then sent to Allframe, which simultaneously photometered all images in a self-consistent manner with regard to centroiding and deblending.

To assemble the final light curves, we first aggregated the Allframe measurements of a given star across all images. Due to the difficulty of obtaining precise flat-field images in space, we cannot calibrate these images using the same methods as would be used for ground-based images. As a result, the light curves of all the stars observed by EPOXI/HRI show variations at the ∼1 % level on a timescale of a few hours, which is the time that it takes for the pointing to drift a distance of the order of the PSF FWHM. Because this is the dominant term in the EPOXI/HRI photometry errors, we bin the data at an interval of 2.4 hr, which seems to remove most, but not all, of the correlations. This gives the light curve shown in the inset in the upper right-hand corner of Figure 1.

Zoom In Zoom Out Reset image size

We had also hoped to get EPOXI/HRI observations after the MOA-2009-BLG-266 microlensing event had returned to its baseline brightness in 2010 March or April. Unfortunately, the EPOXI operations team was busy with preparations for the 2010 November encounter with the comet Hartley 2, so no baseline observations were possible. Therefore, we have determined the baseline brightness in the EPOXI by comparing the EPOXI images to V- and I-band CTIO images taken in 2010 June, when the microlensing event had returned to baseline. The I-band CTIO image is compared to one of the EPOXI frames in Figure 2. Because of the relatively large EPOXI/HRI PSF, we consider only stars in the EPOXI images that have only one counterpart star within a radius of 3'' of the position of the EPOXI star. We also limit our consideration to stars within 0.9 mag of the V − I color of the microlensed target star. There are four stars that satisfy this condition and appear in more than half of the images in which the target star appears. These are the four stars indicated by the green circles in Figure 2. We fit the mean "clear"-filter EPOXI magnitude, C_EPOXI, to the instrumental CTIO magnitudes from the SoDoPHOT reductions, and this yields the following linear relationship between the average EPOXI magnitudes and CTIO magnitudes,

$$\begin{equation} C_{EPOXI} = 0.520 I_{{\rm CTIO}} + 0.480 V_{{\rm CTIO}}. \end{equation} \tag{ 1 }$$

The fit to the magnitude of these four comparison stars gives χ² = 0.22 for 2 degrees of freedom if the uncertainty in the EPOXI magnitudes is assumed to be 0.004 mag. We use this formula to determine the baseline C_EPOXI magnitude, and we add this to the light curve as a final measurement with an assumed uncertainty of 0.01 mag.

Zoom In Zoom Out Reset image size

We note that this calibration procedure for the EPOXI data is probably more reasonably considered to be a procedure to calibrate the source star flux instead of the baseline brightness, which includes the brightness of any unresolved stars blended with the source. But we chose to treat the unmagnified flux estimate as an estimate of the baseline brightness as this is more conservative.

The angular radius of the source star can be determined from its brightness and color, once the effect of interstellar extinction has been removed. We start from the CTIO V- and I-band magnitudes and the IRSF H-band magnitude that have been determined from the best-fit model. The V- and I-band magnitudes have been calibrated to the OGLE-III system (Udalski et al. 2008) and the IRSF H band has been calibrated to the Two Micron All Sky Survey (2MASS; Carpenter 2001).⁷⁴ The comparison between the 2MASS and the IRSF H-band data is subject to complications due to variability and blending, because the 2MASS images, with their 2'' pixels, have significantly worse angular resolution than IRSF. This means that many of the apparent 2MASS "stars" cannot be used for the calibrations because they are actually blended images of two or more stars that are resolved in the IRSF images. This makes calibration of the CTIO H-band images difficult, because of the relatively small 5.5 arcmin² CTIO H FOV. Fortunately, the IRSF FOV is ∼60 arcmin², and it is possible to use over 400 unblended 2MASS stars for the H-band calibration. These calibrations combined with the best-fit light-curve models yield source magnitudes of

$$\begin{eqnarray} H_s =& 13.780 \pm 0.030, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} I_s =& 15.856 \pm 0.030, \end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray} V_s =& 17.677 \pm 0.030, \end{eqnarray} \tag{ 4 }$$

where the uncertainties are almost entirely due to the calibrations (including the uncertainty in the OGLE-III calibration). These magnitudes are indicated by the green dots in the color–magnitude diagrams shown in Figures 3 and 4. The fit uncertainties are ⩽0.005 mag in all three passbands, and the u₀ > 0 model predicts a source that is 0.005 mag brighter than the best-fit u₀ < 0 model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can use the magnitudes from Equations (2) to (4) to determine the source star angular radius, but first we must estimate the foreground extinction. We determine the source radius using the method of Bennett et al. (2010), which is a generalization to three colors of an earlier two-color method (Albrow et al. 2000). Following this procedure, we find the VIH magnitudes of the center of the red clump giant distribution to be

$$\begin{eqnarray} H_{{\rm rc}} =& 13.59 \pm 0.10, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} I_{{\rm rc}} =& 15.73 \pm 0.10, \end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray} V_{{\rm rc}} =& 17.64 \pm 0.10, \end{eqnarray} \tag{ 7 }$$

for stars within 2' of the source star. These are indicated by the red spots in Figures 3 and 4. Assuming a distance to the source of 8.8 kpc (Rattenbury 2007), we can use these red clump magnitudes to estimate the extinction, which we find to be A_H = 0.36, A_I = 1.22, and A_V = 2.14, following an R_V = 2.77 Cardelli et al. (1989) extinction law. These then yield de-reddened magnitudes of H_s0 = 13.42, I_s0 = 14.64, and V_s0 = 15.54. Unlike the case for dwarf stars, the accuracy of the V − I, V − H, and I − H surface-brightness–color relations (Kervella et al. 2004) is similar, but the I − H relation yields an angular radius estimate that is almost completely independent of the reddening law, if it follows the Cardelli et al. extinction law (Cardelli et al. 1989), but this may be due to the fact that the A_I/A_H ratio does not vary much with this extinction law. In any case, all three of these relations imply that the angular radius of the source star is

$$\begin{equation} \theta _\star = 5.2\pm 0.2\ \mu {\rm as}. \end{equation} \tag{ 8 }$$

This and the source radius crossing time of t_⋆ = 0.326 ± 0.007 days imply that the relative lens-source proper motion and Einstein radius are

$$\begin{equation} \mu _{\rm rel} = \theta _\star /t_\star =5.86 \pm 0.26 \ {\rm mas}\ {\rm yr}^{-1}, \end{equation} \tag{ 9 }$$

and

$$\begin{equation} \theta _{\rm E} =\mu _{\rm rel}t_{\rm E} = 0.98\pm 0.04\ {\rm mas}, \end{equation} \tag{ 10 }$$

respectively.

During caustic crossings the lens effectively scans the source star with high angular resolution. As a result, the shape of the light curve reflects the underlying limb darkening of the star (Witt 1995; Bennett & Rhie 1996). Hence, in order to analyze caustic-crossing events such as MOA-2009-BLG-266, one needs to account for the limb darkening appropriately. For the previous planetary microlensing events (Bond et al. 2004; Udalski et al. 2005; Beaulieu et al. 2006; Gould et al. 2006; Gaudi et al. 2008; Bennett et al. 2008; Dong et al. 2009b; Sumi et al. 2010; Janczak et al. 2010; Miyake et al. 2011; Batista et al. 2011), limb darkening has generally been treated within the linear limb-darkening approximation. In some cases, the two-parameter square-root limb-darkening law was used (Dong et al. 2009b; Janczak et al. 2010), even though there has been no indication that the details of the limb-darkening treatment had any noticeable effect on the other model parameters for a planetary microlensing event.

In the case of MOA-2009-BLG-266, there was reason to suspect that the treatment of limb darkening could be important. This is because, as we discuss below in Section 5.1, there are two approximately degenerate microlensing parallax models that have slightly different binary lens parameters. The source crosses the caustics at a slightly different angle for the two models. This suggests that the detailed treatment of limb darkening might have some influence on the difference in χ² between these two degenerate models. As shown by Heyrovsky (2007), using linear limb darkening may introduce photometric errors at the level of 0.01 due to the approximation itself and the choice of method used for computing the linear model coefficients. In order to avoid introducing any such inaccuracies in the analysis of the event, we directly use the limb-darkening profile from a theoretical model atmosphere of the source star, instead of its analytical approximations.

Based on the location of the source star on the color–magnitude diagram, we assume a temperature of T_eff ≈ 4750 K, surface gravity of log g = 2.5, and solar metallicity. We use a model atmosphere from Kurucz's ATLAS9 grid (Kurucz 1996, 1993a, 1993b)⁷⁵ corresponding to these parameters. The raw model data provide values of the specific intensity as a function of wavelength for 17 different positions on the stellar disk. In order to obtain the light-curve-specific limb-darkening profile, we integrate the specific intensity over the relevant filter passband, weighted by the filter transmission, the quantum efficiency of the CCD, and interstellar extinction (see Section 4.1). In order to compute the limb darkening at an arbitrary position on the stellar disk, we interpolate the obtained points using cubic splines with natural boundary conditions (Heyrovsky 2003, 2007). The light-curve modeling code uses pre-computed tabulated intensity values for a sufficiently dense spacing of radial positions on the stellar disk.

This approach avoids a potential source of low-level systematic error without any degrees of freedom to the model. In Table 1 we compare the results of our analysis with those obtained by the usual approach, using linear limb-darkening coefficients from Claret (2000). For the parameters of the source star Claret (2000) provides coefficient values u_λ = 0.7844, 0.7035, 0.6087, 0.4868, 0.4181, and 0.358 for the V, R, I, J, H, and K_s passbands, respectively. Tabulated intensity values give a χ² improvement over the linear approximation of Δχ² = 7.27 for the best-fit static models without orbital motion. Therefore, the tabulated limb-darkening tables fit the data somewhat better, at least for the static lens case, but the implied planetary parameters do not change significantly.

The basic parameters for planetary events like MOA-2009-BLG-266 are straightforward to determine, as a reasonable estimate can be made from the single-lens parameters (found from a fit with the planetary deviation excluded) and inspection of the light curve (Gould & Loeb 1992). The main feature of the deviation is the half-magnitude decrease in magnification centered at HJD' = 5086.5. This indicates that a planet is perturbing the minor (saddle) image created by the stellar lens and that the star–planet separation is less than the Einstein radius. Such a light curve cannot be mimicked by non-planetary perturbations (Gaudi & Gould 1997). The basic planetary parameters can then be estimated following the arguments given in Sumi et al. (2010). In practice, this is not how the parameters were determined, however.

We model the data using standard methods (Bennett 2010; Dong et al. 2006) to extract the precise parameters and uncertainties of the light curve fit. It is convenient to describe microlensing events in terms of the Einstein ring radius, $R_E = \sqrt{(4GM_L/c^2)D_Sx(1-x)}$, which is the radius of ring image seen when the source and (single) lens are in perfect alignment. Here M_L is the lens system mass, x = D_L/D_S, and D_L and D_S are the lens and source distances. Microlensing by a single lens, such as an isolated star, is described by three parameters: the Einstein radius crossing time, t_E, and the time, t₀, and distance (with respect to R_E), u₀, of closest alignment between the source and lens center of mass. Planetary microlensing events require three additional parameters: the planet/star mass ratio, q, the star–planet separation, s, in units of R_E, and angle of the source trajectory with respect to the star–planet axis, θ. The source radius crossing time, t_⋆, is also required because, like most planetary events, MOA-2009-BLG-266 has sharp light curve features that resolve the angular size of the source star.

The MOA and Canopus data for the event were modeled immediately upon the detection of the planetary perturbation using the method of Bennett (2010), supplemented with the addition of the hexadecapole approximation (Pejcha & Heyrovsky 2009; Gould 2008). This found the basic solution that we present here, plus a disfavored alternative s > 1 solution, which was excluded within hours when the planetary deviation data from South Africa, Israel, and Chile became available. The s < 1 solution was refined as more data came in, and the two degenerate solutions that we present here emerged when microlensing parallax was added to the modeling.

We also conducted a blind search of parameter space using the approach of Dong et al. (2006), where the binary parameters s, q, and θ are fixed at a grid of values, while the remaining parameters are allowed to vary so that the model light curve results in minimum χ² at each grid point. A Markov Chain Monte Carlo (MCMC) method was used for χ² minimization. Then, the best-fit model is obtained by comparing the χ² minima of the individual grid points.

The modeling indicates that the perturbation of MOA-2009-BLG-266 is produced by the crossing of a clump giant source star over the planetary caustic produced by a low-mass planet. As we discuss below in Section 5.3, the orbital motion of the planet is not needed to describe the light curve. Assuming a static lens system, the best-fit values of the planet/star mass ratio and normalized star–planet separation are q = 5.425 × 10⁻⁵ and s = 0.91422, respectively. The values of other lensing parameters are listed in Table 2. This table also lists the best-fit parameters for fits including orbital motion. The inclusion of orbital motion improves χ² by Δχ² = 3.24 for 2 fewer degrees of freedom, but it also changes the planet/star mass ratio and normalized star–planet separation to q = 5.815 × 10⁻⁵ and s = 0.91429. As discussed below in Section 5.4, the inclusion of microlensing parallax adds a parameter degeneracy that takes u₀ → −u₀ and θ → −θ, which corresponds to a reflection of the lens plane with respect to the geometry of Earth's orbit. We find that the model with u₀ > 0 is favored by Δχ² = 13.39 for static models and Δχ² = 6.31 for models with orbital motion. The χ² contribution of the individual data sets for the best models with and without orbital motion is shown in Table 1. We note that this mass ratio is the lowest of planets yet to be discovered by the microlensing method.

Table 2. Best-fit Parameters

Parameter	Units	Best Orb.	MCMC Orb.	Best No Orb.	MCMC No Orb.
χ2/dof		6050.17		6053.41
tE	days	61.447	61.47 ± 0.40	61.612	61.54 ± 0.36
t0	HJD'	5093.257	5093.257 ± 0.083	5093.298	5093.285 ± 0.051
u0		0.13158	0.1315 ± 0.0008	0.13129	0.1314 ± 0.0008
s		0.91429	0.91434 ± 0.00036	0.91425	0.91421 ± 0.00036
q	10−5	5.815	5.63 ± 0.25	5.425	5.45 ± 0.07
θ	rad	2.2677	2.2692 ± 0.0034	2.2715	2.2708 ± 0.0035
t⋆	days	0.33140	0.3262 ± 0.0068	0.32152	0.3216 ± 0.0012
πE		0.2094	0.223 ± 0.037	0.2395	0.232 ± 0.038
ϕE	rad	0.760	0.74 ± 0.16	0.640	0.70 ± 0.14
πE, N		0.1517	0.1665 ± 0.0503	0.1921	0.1795 ± 0.0493
πE, E		0.1442	0.1439 ± 0.0035	0.1430	0.1435 ± 0.0035
$\dot{s}_x$	10−3 day−1	−0.34	−0.25 ± 0.15
$\dot{s}_y$	10−3 day−1	−2.05	−0.84 ± 1.49
Torb	days	2380

Note. Parameters are given in an inertial geocentric frame fixed at HJD' = 5086 (HJD' = HJD − 2450, 000).

Download table as:  ASCIITypeset image

Figure 1 shows the best-fit model curve compared to the light curve data. Figure 5 compares the planetary caustic geometry to the source size and trajectory. The two triangular-shaped caustics indicate a minor image caustic perturbation, which is seen when the star–planet separation is less than the Einstein radius (s < 1). The strongest feature in such a minor image caustic crossing event is the large decrease in magnification at HJD ∼ 2, 455, 086.5 when the source is between the caustics and the minor image is essentially destroyed. This feature is surrounded by two positive light curve bumps caused by the source passing over a caustic or passing in front of the cusps. There are no known non-planetary light curve perturbations that can produce such a feature in the light curve (Gaudi & Gould 1997). For MOA-2009-BLG-266, the local light curve minimum between the caustic crossings has a short duration of ∼3.7 hr, which is much smaller than the caustic crossing durations of >20 hr. This is due to the fact that separation between the two triangular minor image caustics is only slightly larger than the diameter of the source star.

Zoom In Zoom Out Reset image size

Like most low-mass planetary microlensing events, MOA-2009-BLG-266 can be modeled well without including any orbital motion of the planet about the host star. But, while we have not measured the light curve precisely enough to measure orbital motion parameters, the planetary orbital motion does influence the precision to which the parameters can be measured from the light curve. In particular, orbital motion allows the planetary caustic to move with respect to the center of mass of the system. Thus, the planetary caustic can be either larger or smaller than the value determined from static lens models, so the mass ratio is not measured as precisely as the static models would imply. In addition, the source radius crossing time is also determined by the time it takes the source to cross the sharp light curve features of the planetary caustic, so this also depends on the orbital motion of the planet and is less precisely determined than the static models would imply.

We include orbital motion using the parameterization used for the analysis of the two-planet event OGLE-2006-BLG-109Lb,c (Bennett et al. 2010), with the x-axis defined by the vector separating the star and planet at HJD = 2, 455, 086. The main orbital motion parameters are $\dot{s}_x$ and $\dot{s}_y$, which describe the instantaneous planet velocity at the time HJD = 2, 455, 086. This parameterization also includes the orbital period, T_orb, but this parameter has a very small effect on the χ² value if it is in the range of physically reasonable values. So, for many of our calculations, we have left T_orb fixed at a physically reasonable value. Independent calculations with a slightly different orbital motion parameterization (Skowron 2011) reached identical conclusions.

The effect of the orbital motion on the other light curve parameters can be seen in Table 2, which shows the parameters and error bars for models with and without orbital motion. The inclusion of orbital motion shifts the values of q and t_⋆ significantly, by 2.6 and 3.8 times the error bars of the static solution, respectively. Orbital motion also increases the error bar on q by a factor of 3.6 and the error bar on t_⋆ by a factor of 5.7, but the error bars on the other parameters do not change significantly, except for the error bar on t₀, which grows by a factor of 1.6. The error bars on the implied physical parameters, shown in Table 3, also do not change very much when orbital motion is included. However, the central values of the physical parameters do change by as much as 0.4σ.

While the light curve has not been measured precisely enough to significantly constrain the orbital motion, the orbital motion can be constrained with the requirement that the planet be bound to its host star. Such a constraint requires that the mass and distance of the host star be known, but the light curve does provide this information as shown below in Equations (13) and (14). The light curve parameters include the transverse host-star–planet separation, s, and the transverse components of the planet velocity, $\dot{s}_x$ and $\dot{s}_y$. One option is to enforce a model constraint that the orbital motion parameters are consistent with a physical circular orbit. This is equivalent to imposing a constraint on the distance to the source (Bennett et al. 2010). We take this limit to be D_S = 8.8 ± 1.2 kpc, which is somewhat larger than the measured scatter of the distance of bulge clump giant stars at the Galactic longitude of this event (Rattenbury 2007) because we do not wish to enforce the circular orbit model too strictly. This constraint has been employed for the best-fit model shown in Table 2. But, such a constraint is inconvenient to use in our MCMC calculations to determine the distribution of allowed light curve and physical parameters, as it makes it more difficult to obtain well-sampled Markov Chains.

A slightly weaker, but more general, constraint can be obtained by requiring that the orbital velocity not be too large to allow the planet to be gravitationally bound to the star, since the probability of lensing by a planet not bound to the lens star is ≲ 10⁻⁸. The transverse velocity components allow us to compute lower limits on the planetary kinetic energy and gravitational binding energy (or an upper limit on the absolute value of the binding energy). The requirement that the total energy < 0 yields an upper limit on the transverse planet velocity (Dong et al. 2009b):

$$\begin{equation} \dot{s}_x^2 + \dot{s}_y^2 \le {2GM_L\over d_\perp R_E^2} = {2GM_L\over s (\theta _E D_L)^3} \, \end{equation} \tag{ 11 }$$

where d_⊥ is the transverse star–planet separation. The R³_E = (θ_ED_L)³ factor in the denominator is needed because the planet–star separation, s, and transverse velocity components use the Einstein radius as their unit of length.

The constraint, Equation (11), has been used in all of our MCMC calculations to determine the allowed parameter distributions, and we have also added a Δχ² = 4 penalty to potential MCMC links with $\dot{s}_x^2 + \dot{s}_y^2$ of more than half the upper limit in Equation (11). Such parameter sets are unlikely because they require a kinetic energy higher than the average value and small values for the separation and velocity along the line of sight. Attempts to find best solutions with the Equation (11) in place of the circular orbit, source distance constraint did not yield better solutions than the one shown in Table 2.

The microlens parallax is defined by the ratio of Earth's orbit to the physical Einstein radius projected on the observer plane, $\tilde{r}_{\rm E}$, i.e.,

$$\begin{equation} \pi _{\rm E}={{\rm AU}\over \tilde{r}_{\rm E}}. \end{equation} \tag{ 12 }$$

Lens parallaxes are usually measured from the deviation of the light curve from those of standard (single or binary) lensing events due to the deviation of the source trajectory from a straight line caused by the orbital motion of the Earth around the Sun (Gould 1992; Alcock et al. 1995). But it is also possible to detect microlensing parallax using observations from a spacecraft in a heliocentric orbit (Refsdal 1966; Dong et al. 2007), and such satellite observations have the potential to significantly increase the number of events for which the microlensing parallax effect may be detected.

For the event MOA-2009-BLG-266, the parallax effect is firmly detected. We find that the χ² difference between the (static) best-fit models with and without the parallax effect is Δχ² = 2789.3, which implies that microlensing parallax is detected at the ∼53σ level. The difference between the parallax and non-parallax models can be seen in Figure 1, where the best-fit model is plotted as a solid black curve and the best-fit non-parallax model is the gray dashed curve. Most of the parallax signal comes from the data outside of the planetary deviation. The light curve without parallax lies below the observed data prior to the planetary perturbation and above the data after the light curve peak. Most of the signal comes from the MOA data, but good coverage of the global light curve shape from CTIO and Canopus has enabled these light curves to contribute to the parallax signal.

There are a number of degeneracies that often affect the parallax parameters of an event. For events with parallax effect, a pair of source trajectories with the impact parameter and source trajectory angle of (u₀, θ) and (− u₀, −θ) can yield degenerate solutions (Smith et al. 2003). Without parallax, this transformation is a trivial redefinition of parameters, but with parallax, we have the reference frame of Earth's orbit, which allows us to distinguish between two solutions that differ by a reflection of the lens plane. For single-lens events with t_E as small as ∼60 days, there is usually an additional degeneracy known as the jerk-parallax degeneracy (Gould 2004), but the additional light curve structure due to the planet removes some of the parameter degeneracy. As a result, the (u₀, θ)↔(− u₀, −θ) degeneracy and the jerk-parallax degeneracy are replaced by a single-parameter degeneracy that makes the (u₀, θ)↔(− u₀, −θ) and changes the parallax parameter, ${\bm {\pi _E}}$. This degeneracy yields the two solutions with similar parameters, but when orbital motion of the planet is ignored, there are significant differences in the other model parameters besides u₀ and θ. With no orbital motion, the u₀ > 0 solution is favored by Δχ² = 13.39. Formally, this is quite significant as the u₀ < 0 solution would be disfavored by a formal probability of $e^{-\Delta \chi ^2/2} \approx 0.0012$, but we must also consider possible systematic errors that may influence the Δχ² value between the u₀ > 0 and u₀ < 0 solutions, as well as the orbital motion of the planet.

This concern about possible systematic errors was the reason for the careful limb-darkening treatment described in Section 4.3. In addition, we also considered a number of different photometric reductions of the data sets that contribute the most to the detection of the parallax signal. These were the MOA, Canopus, and CTIO I- and H-band data sets. With these different photometric reductions, we found that the u₀ > 0 solution was always favored by a similar Δχ² difference, although the SoDoPHOT reduction of the CTIO I-band data set and the DoPHOT reduction of the CTIO H-band data set favored the u₀ < 0 solution by a somewhat larger Δχ² difference. Table 1 indicates that the χ² difference between the u₀ > 0 and u₀ < 0 solutions is spread over a number of different data sets.

The inclusion of the orbital motion parameters discussed in Section 4.3 has a significant effect on the (u₀, θ)↔(− u₀, −θ) degeneracy. When the orbital motion parameters are included, the best-fit (u₀ > 0) model improves its χ² value by Δχ² = 3.24, but the χ² improvement for the u₀ < 0 models is even greater, Δχ² = 10.34, so that the χ² difference between the u₀ < 0 and u₀ > 0 solutions drops to Δχ² = 6.29 when the planetary orbital motion parameters are included in the models. But an even more significant difference is that the differences between the other parameters for these degenerate models largely disappear when orbital motion is included. The added degrees of freedom provided by the orbital motion parameters appear to be larger than the light curve difference enforced by the (u₀, θ)↔(− u₀, −θ) transformation. As a result, once orbital motion is included, the (u₀, θ)↔(− u₀, −θ) degeneracy has no obvious effect on the determination of the physical parameters of the event.

Figure 6 shows the distribution of parallax parameters, (π_{E, N}, π_{E, E}), or equivalently, (π_E, ϕ_E), found by our MCMC simulations. This plot includes 11 separate MCMC chains with a total of 593,000 links as discussed in Section 5.5. The distributions for both solutions are highly elongated along the π_{E, N} axis. This is due to the fact that Earth's acceleration is almost entirely in the east–west direction, when projected on the plane perpendicular to the line of sight to the Galactic bulge. Figure 6 also shows contours of constant π_E, which are labeled by the (approximate) corresponding lens mass. The lens mass depends on the angular Einstein radius, which our MCMC calculations determine to be θ_E = 0.98 ± 0.04 mas. However, these mass contours in Figure 6 are only approximate because they do not include any correlations between the π_E and θ_E values. These correlations are properly incorporated into our MCMC calculations, which yield the host star and planet masses, M_⋆ = 0.56 ± 0.09 M_☉ and m_p  =  10.4 ± 1.7 M_⊕, located at a distance of D_L = 3.04 ± 0.33 kpc. Assuming a random orientation of the orbit, we estimate a semimajor axis of a = 3.2^+1.9 _{− 0.5} AU. If we assume a standard position for the snow line, ∼2.7(M/M_☉) AU (Kennedy & Kenyon 2008), then the planet orbits at about twice the distance of the snow line, similar to the position of Jupiter in our own solar system. Thus, the planet might be considered to be a "failed Jupiter" core as predicted by the core accretion theory (Thommes et al. 2008; Ida & Lin 2005), in which the rock–ice core only reaches ∼10 M_⊕ after the hydrogen and helium gas in the protoplanetary disk has dissipated.

Zoom In Zoom Out Reset image size

In principle, any orbital parallax signal can be mimicked by the so-called xallarap, i.e., orbital motion of the source about a companion (Griest & Hu 1992; Han & Gould 1997). However, this would require very special orbital parameters, basically mimicking those of the Earth (Smith et al. 2002). We search for such xallarap solutions over circular orbits, i.e., with three additional free parameters (orbital phase, inclination, and period). We find a χ² improvement of 3.9 relative to the parallax solution for 3 degrees of freedom, or 3.4 with the period fixed at P = 1 yr (2 degrees of freedom). These improvements have no statistical significance and are to be compared with the improvement of Δχ² = 2789.3 for the parallax solution relative to the no-parallax solution. Therefore, we conclude that the light curve distortions are due to parallax rather than xallarap (Poindexter et al. 2005).

Our final fit parameters are listed in Table 2. The parameters π_{E, N} and π_{E, E} are the north and east components of the microlensing parallax vector, ${\bm {\pi _E}}$. The uncertainty in π_E,N is an order of magnitude larger than the uncertainty in π_{E, E} because the projected acceleration of the Earth is largely in the east–west direction during this event.

Uncertainties in the parameters have been determined by a set of 11 MCMC runs with a total of approximately 593,000 links. Eight of these MCMC runs have been in the vicinity of the u₀ > 0 solution and the other three have been in the vicinity of the u₀ < 0 local χ² minimum. Due to the χ² difference, Δχ² = 6.29 between these local χ² minima, we include a weight factor to our sums over the Markov chains so that the disfavored u₀ < 0 solutions are disfavored by an amount appropriate for their χ² difference, $e^{-\Delta \chi ^2} = 0.043$. The mean parameter values for these solutions and their uncertainties are shown in Table 2. For most parameters, these are given by the weighted averages over the 11 Markov chains, but for u₀ and θ, we have included only the 8 Markov chains with u₀ > 0 and θ > 0. Due to the large difference in these parameters in the vicinity of the two solutions, the mean values would be values that are inconsistent with both solutions if we had used both solutions for these sums. For the remaining parameters, except $\dot{s}_x$ and $\dot{s}_y$, the parameter distributions for the vicinities of the two local χ² minima are nearly identical.

The source radius crossing time, t_⋆, is an important parameter because it helps to determine the angular Einstein radius, θ_E = θ_⋆t_E/t_⋆, as discussed in Section 4.2. When this is combined with the measurement of the microlensing parallax signal, we are able to determine the mass of the lens system (Gould 1992),

$$\begin{equation} M_L = {\theta _E c^2 {\rm AU}\over 4G \pi _E} = {\theta _E \over (8.1439\,{\rm mas})\pi _E} M_\odot \approx 0.57\,M_\odot, \end{equation} \tag{ 13 }$$

if we assume that the favored parameters of the best-fit (u₀ > 0) solution are correct. The lens system distance can also be determined as

$$\begin{equation} D_L={{\rm AU}\over \pi _{\rm E}\theta _{\rm E}+\pi _S} \approx 3.2 \ {\rm kpc}, \end{equation} \tag{ 14 }$$

assuming that the distance to the source, D_S = 1/π_S, is known. Note that these values from the best-fit solution are not identical to the central values from our full MCMC analysis, although they are very close.

In order to determine the physical parameters of this planetary lens system, it is important to include the effects of correlations of the parameters and the uncertainties in external measurements, such as the determination of θ_⋆. Such a calculation is easily done with MCMC simulations. As discussed in Sections 5.4 and 5.5, we have run 11 MCMCs in the vicinity of both the degenerate u₀ > 0 and u₀ < 0. The distribution of the parallax parameters for these solutions is shown in Figure 6. The gray dashed circles in this figure show the contours of constant π_E, which correspond to contours of constant mass by Equation (13). However, this correspondence is only approximate because θ_E is slightly correlated with π_E.

These MCMC simulations can also be used to determine the physical parameters of the host star and its planet, which are summarized in Table 3. This is essentially a Bayesian analysis, but the only non-trivial prior that we impose is the assumption that the orbital orientation is random, which is used to estimate the semimajor axis, a, based on the measured two-dimensional separation in the plane of the sky. If planets are much more common at very small or very large separations, then the planetary detection would imply a bias that violates this assumption. However, the available evidence indicates that planet frequency does not have a sharp dependence on the semimajor axis, so this assumption seems reasonable. Our MCMC calculations assumed a fixed distance of 8.8 kpc to the source star, due to its position on the more distant end of the Galactic bar. We have adjusted the uncertainties in Table 3 to include the 5% spread in the distance to bulge clump stars measured in this direction (Rattenbury 2007; although the effect is quite small.) The probability distributions for the host star and planet masses and distance (M_⋆, m_p, and D_L) are nearly Gaussian, so they are described well simply by their mean values and dispersions. This is not the case for our estimate of the semimajor axis, a, which has a 2σ (95% confidence level) range of 2.3–12.9 AU.

We therefore conclude that the MOA-2009-BLG-266Lb planet is a ∼10 M_⊕ planet at a separation of ∼3 AU. In the core accretion model of planet formation, the snow line is an important location where the density of solid material jumps by about a factor of five due to the condensation of ices (mostly water ice; Ida & Lin 2005; Lecar 2006; Kennedy et al. 2006; Kennedy & Kenyon 2008; Thommes et al. 2008). Assuming a standard position for the snow line, ∼2.7(M/M_☉) AU, we find that the planet is located at about twice the distance of the snow line, similar to the position of Jupiter in our own solar system. It is therefore a prime candidate to be a "failed Jupiter core," which grew by the accumulation of solid material to ∼10 M_⊕, but was unable to grow into a gas giant by the accretion of hydrogen and helium because the protoplanetary disk had lost its gas before the planetary core was massive enough to accrete it efficiently.

The mass measurements of the planet and host star given in Table 3 have uncertainties of about 16%, which is dominated by the uncertainty in π_{E, N}. This uncertainty can be reduced to 5–10 yr, hence when the source and planetary host stars have separated enough to allow their relative proper motion to be measured (Bennett et al. 2007). Since ${\bm {\pi _E}}$ is parallel to the lens–source proper motion, this will reduce the uncertainty in π_E to a value much closer to the 2.4% uncertainty in π_{E, E}, which should reduce the uncertainties in the star and planet masses to <5%. Our existing VLT/NACO observations indicate that the combined H-band flux of the source and host star is H = 13.77 ± 0.05, which is consistent with our prediction that the host star should be ∼75 times fainter than the source and indicate no neighbor stars that might interfere with the detection of the source–host-star relative motion. Therefore, we expect that it will be feasible to improve these mass measurements in the future.

We find a host star of mass M_⋆ = 0.56 ± 0.09 M_☉ orbited by a planet of mass m_p = 10.4 ± 1.7 M_⊕, located at a distance of D_L = 3.04 ± 0.33 kpc. Assuming a random orientation of the orbit, we estimate a semimajor axis of a = 3.2^+1.9 _{− 0.5} AU and an orbital period of P = 7.6^+7.7 _{− 1.4} yr. If we assume a standard position for the snow line, ∼2.7(M/M_☉) AU (Kennedy & Kenyon 2008), as indicated in Figure 7, then the planet orbits at about twice the distance of the snow line, similar to the position of Jupiter in our own solar system. However, the planet's mass of ∼10 M_⊕ is close to the critical mass predicted by core accretion theory (Thommes et al. 2008) where it has exhausted the nearby supply of solid material and begins the slow, quasistatic gas accretion phase. So, MOA-2009-BLG-266Lb fits the theoretical predictions for a large population of "failed gas giant" core (Laughlin et al. 2004) planets, which would have had their growth terminated by the loss of gas from the protoplanetary disk prior to the rapid gas accretion phase. Indeed, the distribution of planets found by microlensing, shown in Figure 7, seems to confirm this prediction, as the detection efficiency corrected planetary mass function rises steeply, as ∼q^{−0.7 ± 0.2} toward lower mass ratios (Sumi et al. 2010). However, a much sharper comparison to theory can be made with mass measurements of these planets and their host stars. Some theoretical treatments suggest a relatively sharp feature in the mass function at ∼10 M_⊕, and the low-mass end of the exoplanet mass function is likely to depend on the host star mass. MOA-2007-BLG-192Lb (Kubas et al. 2011) is the only other cold, low-mass planet with a host star mass measurement, but the planet mass is weakly constrained due to a poorly sampled light curve.

Zoom In Zoom Out Reset image size