A Giant Planet beyond the Snow Line in Microlensing Event OGLE-2011-BLG-0251

We present the analysis of the gravitational microlensing event OGLE-2011-BLG-0251. This anomalous event was observed by several survey and follow-up collaborations conducting microlensing observations towards the Galactic Bulge. Based on detailed modelling of the observed light curve, we find that the lens is composed of two masses with a mass ratio q=1.9 x 10^-3. Thanks to our detection of higher-order effects on the light curve due to the Earth's orbital motion and the finite size of source, we are able to measure the mass and distance to the lens unambiguously. We find that the lens is made up of a planet of mass 0.53 +- 0.21,M_Jup orbiting an M dwarf host star with a mass of 0.26 +- 0.11 M_Sun. The planetary system is located at a distance of 2.57 +- 0.61 kpc towards the Galactic Centre. The projected separation of the planet from its host star is d=1.408 +- 0.019, in units of the Einstein radius, which corresponds to 2.72 +- 0.75 AU in physical units. We also identified a competitive model with similar planet and host star masses, but with a smaller orbital radius of 1.50 +- 0.50 AU. The planet is therefore located beyond the snow line of its host star, which we estimate to be around 1-1.5 AU.


Introduction
Gravitational microlensing is one of the methods that allows us to probe the populations of extrasolar planets in the Milky Way, and has now led to the discoveries of 16 planets 1 , several of which could not have been detected with other techniques (e.g. Beaulieu et al. 2006, Gaudi et al. 2008, Muraki et al. 2011. In particular, microlensing events can reveal cool, low-mass planets that are difficult to detect with other methods. Although this method presents several observational and technical challenges, it has recently led to several significant scientific results. Sumi et al. (2011) analysed short time-scale microlensing events and concluded that these events were produced by a population of Jupiter-mass free-floating planets, and were able to estimate the number of such objects in the Milky Way. Cassan et al. (2012) used 6 years of observational data from the PLANET collaboration to build on the work of Gould et al. (2010) and Sumi et al. (2011), and derived a cool planet mass function, suggesting that, on average, the number of planets per star is expected to be more than 1.
Modelling gravitational microlensing events has been and remains a significant challenge, due to a complex parameter space and computationally demanding calculations. Recent developments in modelling methods (e.g. Cassan 2008;Kains et al. 2009Kains et al. , 2012Bennett 2010;Ryu et al. 2010;Bozza et al. 2012), however, have allowed microlensing observing campaigns to optimise their strategies and scientific output, thanks to real-time modelling providing prompt feedback to observers as to the possible nature of ongoing events.
In this paper we present an analysis of microlensing event OGLE-2011-BLG-0251, an anomalous event discovered during the 2011 season by the OGLE collaboration and observed intensively by follow-up teams. In Sec. 2, we briefly summarise the basics of relevant microlensing formalism, while we discuss our data and reduction in Sec. 3. Our modelling approach and results are outlined in Sec. 4; we translate this into physical parameters of the lens system in Sec. 5 and discuss the properties of the planetary system we infer.

Microlensing formalism
Microlensing can be observed when a source becomes sufficiently aligned with a lens along the line of sight that the deflection of the source light by the lens is significant. A characteristic separation at which this occurs is the Einstein ring radius. When a single point source approaches a single point lens of mass M with a projected source-lens separation u, the source brightness is magnified following a symmetric "point source-point lens" (PSPL) pattern which can be parameterised with an impact parameter u 0 and a timescale t E , both expressed in units of the angular Einstein radius (Einstein 1936), where G is the gravitational constant, c is the speed of light, and D S and D L are the distances to the source and the lens, respectively, from the observer. The timescale is then t E = θ E /µ, where µ is the lens-source relative proper motion. Therefore the observable t E is a degenerate function of M, D L and the source's transverse velocity v ⊥ , assuming that D S is known. However, measuring certain second-order effects in microlensing light curves such as the parallax due to the Earth's orbit allows us to break this degeneracy and therefore measure the properties of the lensing system directly. When the lens is made up of two components, the magnification pattern can follow many different morphologies, because of singularities in the lens equation. These lead to source positions, along closed caustic curves, where the lensing magnification is formally infinite for point sources, although the finite size of sources means that, in practice, the magnification gradient is large rather than infinite. A point-source binary-lens (PSBL) light curve is often described by 6 parameters: the time at which the source passes closest to the center of mass of the binary lens, t 0 , the Einstein radius crossing time, t E , the minimum impact parameter u 0 , which are also used to describe PSPL light curves, as well as the source's trajectory angle α with respect to the lens components, the separation between the two mass components, d, and their mass ratio q. Finite source size effects can be parameterised in a number of ways, usually by defining the angular size of the source ρ * in units of θ E : where θ * is the angular size of the source in standard units.

Observational data
The microlensing event OGLE-2011-BLG-0251 was discovered by the OGLE (Optical Gravitational Lens Experiment) collaboration's Early Warning System (Udalski 2003) as part of the release of the first 431 microlensing alerts following the OGLE-IV upgrade. The source of the event has equatorial coordinates α = 17 h 38 m 14.18 s and δ = −27 • 08 ′ 10.1 ′′ (J2000.0), or Galactic coordinates of (l, b)=(0.670 • , 2.334 • ). Anomalous behaviour was first detected and alerted on August 9, 2011 (HJD∼2455782.5) thanks to real-time modelling efforts by various follow-up teams that were observing the event, but by that time a significant part of the anomaly had already passed, with sub-optimal coverage due to unfavourable weather conditions. The anomaly appears as a two-day feature spanning HJD = 2455779.5 to 2455781.5, just before the time of closest approach t 0 . Despite difficult weather and moonlight conditions, the anomaly was securely covered by data from five follow-up telescopes in Brazil (µFUN Pico dos Dias), Chile (MiNDSTEp Danish 1.54m) New Zealand (µFUN Vintage Lane, and MOA Mt. John B&C), and the Canary Islands (RoboNet Liverpool Telescope). The descending part of the light curve also suffered from the bright Moon, with the source ∼ 5 degrees from the Moon at ∼ 85% of full illumination, leading to high background counts in images and more scatter in the reduced data. We opted not to include data from Mt. Canopus 1m telescope in the modelling because of technical issues at the telescope affecting the reliability of the images, and also excluded the I-band data from CTIO because they also suffer from large scatter, probably due to the proximity of the bright full Moon to the source.
The data set amounts to 3738 images from 13 sites, from the OGLE survey team, the MiNDSTEp consortium, the RoboNet team, as well as the µFUN, PLANET and MOA collaborations in the I, V and R bands, as well as some unfiltered data; data sets are summarised in Table 1 and the light curve is shown in Fig. 1. We reduced all data using the difference imaging pipeline DanDIA (Bramich 2008;Bramich et al. 2013), except for the OGLE data, which was reduced by the OGLE team with their optimised offline pipeline.
For each data set, we applied an error bar rescaling factors a and b to normalise error bars with respect to our best-fit model (see Sec. 4), using the simple scaling relation σ ′ i = a σ 2 i + b 2 where σ ′ i is the rescaled error bar of the i th data point and σ i is the original error bar. The error bar rescaling factors for each data set is given in Table 1. We did not exclude outliers from our data sets, unless we had reasons to believe that an outlier had its origin in a bad observation, or in issues with the data reduction pipeline.

Modelling
We modelled the light curve of the event using a Markov Chain Monte Carlo (MCMC) algorithm with adaptive step size. We first used the "standard" PSBL parameterisation in our modelling, whereby a binary-lens light curve can be described by 6 parameters: those given in Sec. 2, ignoring the second-order ρ * parameter described in that section. For all models and configurations we searched the parameter space for solutions with both a positive and a negative impact parameter u 0 .
We started without including second-order effects of the source having a finite size or parallax due to the orbital motion of Earth around the Sun, and then added these separately in subsequent modelling runs by fitting the source size parameter ρ * , as defined in Sec. 2, and the parallax parameters described below. Both effects led to a large decrease in the χ 2 statistic of the model (> 1000), which could not be explained only by the extra number of parameters.
For the finite-source effect, we additionally considered the limb-darkening variation of the source star surface brightness by modelling the surface-brightness profile as where I 0,ψ is the brightness at the centre of the source, and ψ is the angle between a normal to the surface and the line of sight. We adopt the limb-darkening coefficients based on the source type determined from the dereddened colour and brightness (see Sec. 5.1). The values of the adopted coefficients are c V = 0.073, c I = 0.624, c R = 0.542, based on the catalogue of Claret (2000). Finally, in a third round of modelling, we included both the effects of parallax and finite source size ("ESBL + parallax"). Including these effects together led to a significant improvement of the fit, with ∆χ 2 > 500 compared to the fits in which those effects were added separately. Computing the f -statistic (see e.g. Lupton 1993) for this difference tells us that the probability of this difference occurring solely due to the number of degrees of freedom decreasing by 1 or 2 is highly unlikely. Our best-fit ESBL + parallax model is shown in Fig. 1.
To model the effect of parallax, we used the geocentric formalism (Dominik 1998;An et al. 2002;Gould 2004), which has the advantage of allowing us to obtain a good estimate of t 0 , t E and u 0 from a fit that does not include parallax. This formalism adds a further 2 parallax parameters, π E,E and π E,N , the components of the lens parallax vector π E projected on the sky along the east and north equatorial coordinates, respectively. The amplitude of π E is then Measuring π E in addition to the source size allows us to break the degeneracy between the mass, distance and transverse velocity of the lens system that is seen in Eq. (1). This is because π E also relates to the lens and source parallaxes π L and π S as Using this in Eq.
(1) allows us to solve for the mass of the lens.
As an additional second-order effect, we also consider the orbital motion of the binary lens. Under the approximation that the change rates of the binary separation and the rotation of the binary axis are uniform during the event, the orbital effect is taken into consideration with 2 additional parameters ofḋ andα, which represent the rate of change of the binary separation and the source trajectory angle with respect to the binary axis, respectively. It is found that the improvement of fits by the orbital effect is negligible and thus our best-fit model is based on a static binary lens.
Below we outline our modelling efforts that resulted in fits that were not competitive with our best-fit ESBL + parallax models, and which we therefore excluded in our light curve interpretation.

Xallarap
We attempted to model the effects of so-called xallarap, orbital motion of the source if it has companion (Griest & Hu 1992). Modelling this requires five additional parameters: the components of the xallarap vector, ξ E,N and ξ E,E , the orbital period P , inclination i and the phase angle ψ of the source orbital motion. By definition, the magnitude of the xallarap vector is the semi-major axis of the source's orbital motion with respect to the centre of mass, a S , normalised by the projected Einstein radius onto the source plane,r E = D S θ E , i.e.
The value of a S is then related to the semi-major axis of the binary by where M 1 and M 2 are the masses of the source components.  Fig. 2, we show χ 2 of the fit plotted as a function of the orbital period of the source star. We compare this to the χ 2 statistic of the best parallax fit. We find that xallarap models provide fits competitive with the parallax planetary models for orbital periods P > 1 year. However, the solutions in this range cannot meet the constraint provided by the source brightness. Combining Equations (6) and (7) with Kepler's third law, P 2 = a 3 /(M 1 + M 2 ) yields (Dong et al. 2009) Rearranging this equation for M 2 , and using the fact that M 2/(M 1 + M 2 ) < 1, we can derive an upper limit for the mass of M 2 , In the lower panel of Fig. 2, we show the minimum mass of the source companion as a function of orbital period. The blending constraint means that the source companion cannot be arbitrarily massive, and we use a conservative upper limit for its mass of 3 M ⊙ . With this constraint, we find that xallarap models are not competitive with parallax planetary models, and we therefore exclude the xallarap interpretation of the light curve.

Binary source
We also attempted to model this event as a binary source -point lens (BSPL) event. For this we introduced three additional parameters: the impact parameter of the secondary source component, u 0,2 , and its time of closest approach, t 0,2 , as well the flux ratio between the source components. We note that parallax is also considered in our binary source modelling, for fair comparison to other models. We find that the best binary-source model provides a poorer fit, with χ 2 = 3809, which gives ∆χ 2 ∼ 180 compared to our best planetary model (including parallax, see model D in the following section). Residuals for this model, as well as all other models discussed in this section are shown in Fig. 3 Table 1. Data sets for OGLE-2011-BLG-0251, with the number of data points for each telescope/ filter combination. The rescaling coefficients a and b are also given, with error bars rescaled as σ ′ = a √ σ 2 + b 2 , where σ ′ is the rescaled error bar and σ is the original error bar.

Best-fit models
We searched the parameter space using an MCMC algorithm as well as a grid of (d, q, α) to locate good starting points for the algorithm (see e.g. Kains et al. 2009), over the range −4 < log q < 0 and −1.0 < log d < 2. This encompasses both planetary and binary companions that might cause the central perturbation. In Fig. 4 we present the χ 2 distribution in the d, q plane. We find four local solutions, all of which have a mass ratio corresponding to a planetary companion. We designate them as A, B, C and D; the degeneracy among these local solutions is rather severe, as can be seen from the residuals shown in Fig. 3.
For the identified local minima, we then further refine the lensing parameters by conducting additional modelling, considering higher-order effects of the finite source size and the Earths orbital motion. It is found that the higher-order effects are clearly detected with ∆χ 2 > 500. Best-fit parameter for each of the local minima are given in Table 2, while Fig. 5 shows the geometry of the source trajectories with respect to the caustics for all four minima. We note that the pairs of solutions A and D, and B and C, are degenerate under the well-known d ↔ d −1 degeneracy (Griest & Safizadeh 1998;Dominik 1999); this is caused by the symmetry of the lens mapping between binaries with d and d −1 . Comparing the pairs of solutions, we find that the A-D pair is favoured, with ∆χ 2 > 40 compared to the B-C pair. On the other hand, the degeneracy between the A and D solutions is very severe, with only ∆χ 2 ∼ 7. In Fig. 6, we also show parameter-parameter correlations plots for model D, showing also the uncertainties in the measured lensing parameters.

Lens Properties
In this section we determine the properties of the lens system, using our best-fit model parameters, i.e. our wideconfiguration ESBL + parallax model. We also calculated the lens properties for the competitive close-configuration model, with both sets of parameter values listed in Table 3.

Source star and Einstein radius
We determined the Einstein radius by first calculating the angular size of the source. This can be done by using the magnitude and colour of the source (e.g. Yoo et al. 2004), and empirical relations between these quantities and the angular source size. We start by using the location of the red giant clump (hereafter RC) on our colour-magnitude diagram (Fig. 7) to estimate the reddening and extinction along the line of sight. We use an I-band absolute magnitude for the RC of M I,RC,0 = −0.12 ± 0.09 (Nataf et al. 2012), as well as a colour (V − I) RC,0 = 1.06 ± 0.12 Figure 3. Residual of data, with 1-σ error bars, for the various models considered.  Table 2. Best-fit model parameters and 1-σ error bars for the four identified best binary-lens models including the effects of the orbital motion of the Earth (parallax). MHJD≡HJD-2450000. a for the OGLE data set (Bensby et al. 2011). We compare these values to those on our colour-magnitude diagram (CMD), which we generated using OGLE I− and V − band photometry. From Fig. 7, the location of the RC on our CMD is (I, V − I) RC = (17.19 ± 0.05, 3.45 ± 0.05) so, using a distance modulus of µ = 14.52 ± 0.09, i.e. a distance to the Galactic bulge of 8.0 ± 0.3 kpc (Yelda et al. 2011), we find A I = 2.79 ± 0.10 and E(V − I) = 2.39 ± 0.13.
Using these values, the best-fit value for the magnitude of the source I S = 15.97±0.01, a source colour (V −I) S,0 = 1.15, and the empirical relations of Kervella & Fouqué (2008), we find an angular source radius θ * = 10.41 ± 1.18 µas, or a source star radius of R * = 10.53 ± 1.19 R . This, together with the best-fit value of the source size parameter ρ * , allows us to calculate the size of the Einstein radius, θ E = θ * /ρ * . Using the relevant parameter values, we find θ E = 0.749 ± 0.283 mas. This in turn allows us to calculate the source-lens relative proper motion, µ rel = θ E /t E = 4.28 ± 1.62 mas/yr.

Masses of the Lens Components
Combining Eq. (1) and Eq. (5) allows us to derive an expression for the mass as a function of the parallax vector magnitude π E (defined by Eq. 4): Figure 4. χ 2 map in the d, q plane, showing the location of the four local minima identified by our modelling runs. Out of these, local minima A and D are competitive, with local minima B and C having ∆χ 2 ∼ 50 and 70 respectively, for the same number of parameters. Minima A and D correspond to the close and wide ESBL + parallax models discussed in the text. Different colours correspond to ∆χ 2 < 25 (red), 100 (yellow), 225 (green), and 400 (blue); we note that the χ 2 map is based on the original data, before error-bar normalisation, and therefore the ∆χ 2 levels are slightly different from those given in Table 2. The top panel shows the breadth of our parameter space exploration, encompassing planetary and non-planteray mass-ratio regimes, while the bottom panel shows a zoom on the region where our local minima are located.
Using values found in the previous section, and our bestfit parallax parameter value π E = 0.35 ± 0.05 yields a total lens mass M L = 0.26 ± 0.10 M ⊙ . Using the best-fit mass ratio parameter value of q = (1.92 ± 0.12) × 10 −3 yields component masses of 0.26 ± 0.11 M ⊙ and 0.53 ± 0.21 M J , where M J is the mass of Jupiter.

Distance to the Lens
We can also rearrange Eq. (1) to derive an expression for the distance to the lens D L , Using our parameter values as well as the lens mass derived thanks to our parallax measurement, we find a distance to the lens of D L = 2.57 ± 0.61 kpc. This distance allows us to carry out a sanity check of the lens mass we derived in the previous section. By assuming that the contribution from the blended light comes from the lens, we can derive an upper limit to the I−band lens magnitude M I using our best-fit blending parameter: Figure 5. Source trajectory geometry with respect to the caustics for all four local minima identified in Fig. 4; the source size is marked as a red circle.
where m I,b is the apparent I−band magnitude of the blend, A I,L is the extinction between the observer and the lens, and D L is in kpc. In practice, A I,L ≤ A I since the lens is in front of the source, so we use the extreme scenario where A I,L = A I to derive an upper brightness limit (lower limit on the magnitude) for the lens. We find this to be M I,L = 2.19 ± 0.53 mag, which corresponds to a maximum mass of the lens of M L,max = 1.65 ± 0.23M ⊙ , assuming a main sequence star mass-luminosity relation, and assuming that the secondary lens component (i.e. the planet) does not contribute to the blended light. This is much larger than the value we derived in Sec. 5.2 for the mass of the primary lens component, which suggests that some blending comes from stars near the source rather than from the lens, although it is difficult to quantify this without an estimate of A I,L . Finally, we can also use the distance to the lens and the size of the Einstein ring radius to calculate the projected separation r ⊥ between the lens components in AU. Using our best-fit projected angular separation d = 1.408 ± 0.019, we find a projected (i.e. minimum) orbital radius r ⊥ = 2.72 ± 0.75 AU.
We can compare this to an estimate of the location of the "snow line", which is the location at which water sublimated in the midplane of the protoplanetary disk, i.e. the distance at which the midplane had a temperature of T mid = 170 K (although other studies have noted that this temperature varies; see e.g. Podolak & Zucker 2004). The core accretion model of planet formation predicts that giant planets form much more easily beyond the snow line, thanks to easier condensation of icy material and therefore easier formation of large solid cores in the early stages of the circumstellar disk's evolution. Kennedy & Kenyon (2008) modelled the evolution of the snow line's location, taking into account heating of the disk via accretion, as well as the influence of pre-main sequence stellar evolution. Using a rough extrapolation of their results, we estimate that the snow line (at t = 1Myr) for the planetary host star in OGLE-2011-BLG-0251 is located at around ∼ 1−1.5 AU. We therefore conclude that OGLE-2011-BLG-0251Lb is a giant planet located beyond the snow line, with both of our competitive best-fit models yielding projected orbital radii larger than 1.5 AU.  Table 3. Lens properties derived as detailed in Sec. 5, for both competitive parallax models.
We list all the lens properties in Table 3, both for the best-fit model parameters that we have used above, and for the close-configuration model, for comparison. Lens properties derived using the close-configuration model are very similar to those we found using the wide-configuration model, the only major difference being in the orbital radius. For the close model, we find an orbital radius of 1.50 ± 0.50 AU, which is close to the location of the snow line.

Conclusions
Our coverage and analysis of OGLE-2011-BLG-0251 has allowed us to locate and constrain a best-fit binary-lens model corresponding to an M star being orbited by a giant planet. This was possible through a broad exploration of the parameters both in real time, thanks to the recent developments in microlensing modelling algorithms, and after the source had returned to its baseline magnitude. Various second-order effects, as well as other possible, nonplanetary, interpretations for the anomaly were considered during the modelling process. Based on the best-fit solution, we were able to constrain the masses and separation of the lens components, as well as various other characteristics, thanks to a strong detection of parallax effects due to the Earth's orbit around the Sun, in conjunction with the detection of finite source size effects. We found a planet of mass 0.53 ± 0.21 M J orbiting a lens of 0.26 ± 0.11M ⊙ at a projected radius r ⊥ = 2.72 ± 0.75 AU; the whole Figure 6. Parameter-parameter correlations for our 9 fitted parameters. Colours indicate the limits of the 1, 2, 3, 4 and 5-σ confidence limits for each pairwise distribution. A closer view of the correlation between parallax parameters is shown on the top right inset.
system is located at a distance of 2.57 ± 0.61 kpc. Our competitive second-best model leads to similar properties, but a smaller projected orbital radius r ⊥ = 1.50 ± 0.50. The two best-fit models are competitive and therefore we cannot make a strong claim about which orbital radius is favoured. However, comparing both values of the projected orbital radius to the approximate location of the snow line for a typical star of the mass of the primary lens component, we conclude that OGLE-2011-BLG-0251Lb is a giant planet located around or beyond the snow line. This is in line with predictions from the core accretion model of planet formation, from which we expect large planets to be more numerous beyond the snow line; this is also where microlensing detection sensitivity is at its highest, enabling us to probe a region of planetary parameter space that is difficult to reach for other methods. The location of the total source + blend is indicated by a green asterisk, while the location of the deblended source is marked by a blue filled circle, and that of the blend by a red cross. The dashed lines cross at the location of the Red Clump.
knowledge support from NSF AST-1103471. B.S. Gaudi, A. Gould, and R.W. Pogge acknowledge support from NASA grant NNX12AB99G. The MOA experiment was supported by grants JSPS22403003 and JSPS23340064. TS was supported by the grant JSPS23340044. Y. Muraki acknowledges support from JSPS grants JSPS23540339 and JSPS19340058.