OGLE-2014-BLG-0221Lb: A Jupiter Mass Ratio Companion Orbiting Either a Late-type Star or a Stellar Remnant

We present the analysis of the microlensing event OGLE-2014-BLG-0221, a planetary candidate event discovered in 2014. The photometric light curve is best described by a binary-lens single-source model. Our light-curve modeling finds two degenerate models, with event timescales of t E ∼ 70 days and ∼110 days. These timescales are relatively long, indicating that the discovered system would possess a substantial mass. The two models are similar in their planetary parameters with a Jupiter mass ratio of q ∼ 10−3 and a separation of s ∼ 1.1. Bayesian inference is used to estimate the physical parameters of the lens, revealing that the shorter timescale model predicts 65% and 25% probabilities of a late-type star and white dwarf host, respectively, while the longer timescale model favors a black hole host with a probability ranging from 60% to 95%, under the assumption that stars and stellar remnants have equal probabilities of hosting companions with planetary mass ratios. If the lens is a remnant, this would be the second planet found by microlensing around a stellar remnant. The current separation between the source and lens stars is 41–139 mas depending on the models. This indicates the event is now ready for high-angular-resolution follow-up observations to rule out either of the models. If precise astrometric measurements are conducted in multiple bands, the centroid shift due to the color difference between the source and lens would be detected in the luminous lens scenario.


INTRODUCTION
Gravitational microlensing is a well-known method for the discovery of exoplanets.Since the first discovery of a planet via microlensing in 2004 (Bond et al. 2004;Bennett et al. 2006), the number of exoplanets detected by microlensing has grown to 200 (Akeson et al. 2013).The most likely lens star in a microlensing event is a low mass late type star as these are the most prolific in our Galaxy.A large fraction of planets discovered by microlensing are gas giants 1 .Such planets are less likely to be formed around low mass stars according to the core accretion model of planet formation (Ida & Lin 2005;Burn et al. 2021).Suzuki et al. (2016) gave a statistical analysis of the set of microlensing planets and found the planet-star mass ratio function was best characterized by a broken power-law with the break at q = 1.7 × 10 −4 which corresponds to the mass of Neptune.Comparison of the result to population synthesis models of the core accretion admitted excess in the microlensing planet beyond q ∼ 10 −4 (Suzuki et al. 2018), including a gas giant mass ratio regime.
Microlensing surveys were originally initiated to search for dark matter in the form of MAssive Compact Halo Objects (MACHOs; Paczynski 1986) in which stellar remnants, substellar objects and planets are the expected populations.Identifying isolated black holes as a part of a population of MACHOs was not successful until recently owing to the degeneracy in mass and distance present in the interpretation of most single lens events.The first definitive detection of an isolated black hole was made possible through observations of the astrometric shift of the apparent source position using high resolution imaging of the Hubble Space Telescope (Sahu et al. 2022;Lam et al. 2022;Mróz et al. 2022).
The first detection of a stellar remnant by microlensing was made in 2021 (Blackman et al. 2021).Following the actual discovery in 2010 (Bachelet et al. 2012), follow-up observations with high resolution imaging using the Keck-II telescope were conducted to resolve the source and lens and resulted in no detection of a lumi-1 110 discovered microlensing planets have estimated mass larger than Saturn (Akeson et al. 2013).Note that microlensing has higher detection efficiency for higher mass ratio planets.
nous lens.The lens system was concluded to be composed of a white dwarf and accompanying Jovian planet.So far, these are the only examples of stellar remnants found by microlensing while there are a few more candidates to be confirmed (e.g.Miyake et al. 2012;Shvartzvald et al. 2015).An ambitious program with the Nancy Grace Roman Space Telescope, the NASA's next flagship mission (Spergel et al. 2015), has been proposed to detect isolated black holes, utilizing Roman's high precision photometry and astrometry (Lam et al. 2023).The survey strategy required for detection of isolated black holes is mostly satisfied by the notional design of the Galactic Bulge Time Domain Survey optimized for detection of exoplanets toward the Galactic bulge.Lam et al. (2023) predicts more than 300 isolated black hole detections among which 270 black holes can be characterized by additional daily cadence observations and an astrometric precision of 0.1 mas.The same strategy should also enable the discovery of black hole binaries.These are less common but more likely to be characterized with additional constraints on the mass distance degeneracy from accurate measurement of a projected source size.A future black hole survey using the powerful capability of Roman seems promising; however, we need to wait until its expected launch in 2026.A discussion now of the discoveries of black holes and other compact objects is important for establishing the future prospects of such a survey.
This paper presents analysis of microlensing event OGLE-2014-BLG-0221.Section 2 describes how the event was discovered and observed.Section 3 gives the details of our light curve analysis.In section 4, we investigate the source property using a color magnitude diagram for field stars around the event coordinate.The lens properties are estimated in section 5.Then, we discuss the results in section 6.

OBSERVATION
Microlensing event OGLE-2014-BLG-0221 was detected by the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 2015) Collaboration on 2014 March 6 (HJD ′ ≡ HJD − 2450000 ∼ 6723) as a part of the OGLE-IV survey.The OGLE group conducted the observations with their 1.3m Warsaw tele-scope at the Las Campanas Observatory in Chile.The event coordinate was first reported at (RA, Dec) J2000 = (18h01m12.90s,−27 • 25 ′ 37.2 ′′ ), corresponding to Galactic coordinates (l, b) = (3.044660,−2.195338),where a known faint star is located; however later, a centroid shift of the lensed source was observed while the source was magnified.The coordinate was corrected to (RA, Dec) J2000 = (18h01m12.90s,−27 • 25 ′ 36.35 ′′ ), corresponding to Galactic coordinates (l, b) = (3.044658,−2.195316), and the event was re-identified by the New Objects in the OGLE Sky (NOOS) system (Udalski 2003) as the lensing of a previously unknown object and re-designated as OGLE-2014-BLG-0284.OGLE-2014-BLG-0221 was located in the OGLE field BLG511 and was detected using a nominal cadence of every 60 minutes for their I -band observations.
The event was also independently discovered by the Microlensing Observations in Astrophysics (MOA; Bond et al. 2001;Sumi et al. 2003) Collaboration as a part of the MOA-II survey on 2014 March 9 and designated as MOA-2014-BLG-069.The event was located in MOA field gb10.This field was one which was observed in the survey with a cadence of 15 minutes.The MOA group uses their 1.8m MOA-II telescope and 2.2 deg 2 wide field of view camera, MOA-cam3, at the Mount John Observatory in New Zealand (Sako et al. 2008).
The OGLE photometric data were obtained mostly in the standard Kron-Cousins I -band and occasionally in the standard Johnson V -band in order to extract color information of the source star.The MOA photometric data were obtained in the designated MOA-Red band, equivalent to a combined band of the standard Cousins R and I. MOA also occasionally observes in JohnsonVband; however unfortunately, no data were taken in 2014.OGLE and MOA reduced the data with their own pipelines (Udalski 2003;Bond et al. 2001) based on their implementation of the difference image analysis method (Alard & Lupton 1998).
The pipelines often underestimate uncertainties of photometric data for a stellar dense region.For this reason, we rescaled the errors to account for low level unknown systematics so that the reduced χ 2 (or χ 2 /d.o.f.) equals 1, following the standard procedure detailed in (Bennett et al. 2008;Yee et al. 2012).The rescaling formula is where σ ′ i is the rescaled error, σ i is the error before rescaling, and k and e 2 min are the rescaling parameters.Data sets and rescaling parameters adopted in the analysis are listed in Table 1.We used the optimized photometry of OGLE-2014-BLG-0221 for OGLE data.The number of data points. b The rescaling parameters.

LIGHT CURVE ANALYSIS
The light curve of OGLE-2014-BLG-0221 in Figure 1 shows a clear anomaly from HJD ′ ∼ 6738 to 6746, a deviation from the symmetric Paczynski curve (Paczynski 1986) with a caustic crossing like shape, indicating the presence of a binary lens system.We modeled the light curve with the following formularization.
The timescale of the event is defined as where µ rel is the relative source-lens proper motion and is the angular Einstein radius as a function of the constant κ = 4G/c 2 , the total lens mass M , and the relative source-lens parallax The time at the source-lens closest approach and the impact parameter are parameterized as t 0 and u 0 .The ratio of the angular source radius θ * and the angular Einstein radius is also measured as For the case of a binary lens event, the mass ratio of the primary and companion q, the projected angular separation of them in units of the angular Einstein radius s, and the sky projected source incident angle α relative to the binary lens axis are included as additional parameters.We considered annual parallax effect and linear lens orbital motion modeled by two additional parameters for each, (π E,N , π E,E ) and (γ ∥ = ds dt , γ ⊥ = dα dt ), where π E,N and π E,E are the north and east components of the microlensing parallax vector (Gould 2004) and ds dt and dα dt are the projected linear and angular motions of the lens companion, expressed as the rate of change in s and α, around the primary (Skowron et al. 2011).
As boundary conditions of a gravitationally bound lens companion, we constrain the orbital motion parameters such that the ratio of the projected kinetic to potential energy of the lens ( KE PE ) ⊥ does not go beyond 0.5, under an assumption of a circular orbit (Dong et al. 2009).Source orbital motion, another higher order effect called xallarap (Poindexter et al. 2005), arising if the source star has a companion, is known to mimic the effects of annual parallax and lens orbital motion.While we also modeled this effect with seven additional parameters (ξ E,N , ξ E,E , R.A. ξ , decl.ξ , P ξ , e ξ , T peri ), we will not go into detail in this paper since we found no solid evidence for xallarap.We applied a linear limb darkening model to the source star in our modeling.With T eff ∼ 5500 K estimated from the source color (see Section 4) and an assumption of surface gravity log g = 4.5 cm/s 2 and solar metallicity log[M/H] = 0, the limb darkening coefficients u I = 0.5189, u V = 0.6854 and u MOA−Red = 0.6052 were determined from the tables of Claret & Bloemen (2011).
We searched for the best-fit binary lens model by exploring the parameter space using our modeling software (Sumi et al. 2010) which is based on the Markov Chain    Monte Carlo method (Verde et al. 2003) and the image centered ray-shooting method (Bennett & Rhie 1996).
In order to explore a wide range of parameter values, we start from a grid search in which we divide the (q, s, α) parameter space into 11×22×40 uniformly spaced grids with logq ∈ [−4, 0], logs ∈ [−0.5, 0.55], α ∈ [0, 2π).We search for a best-fit model at each grid by fixing (q, s, α) and allowing other parameters free.Among the 9680 models, 100 models with the lowest χ 2 values are then refined to explore the global minimum by releasing the fixed (q, s, α) to vary.We finally exclude all the models that exceed the threshold of ∆χ 2 = 100 from the best-fit lowest χ 2 model, through which we find two degenerate models remain.
Figure 2 shows the two degenerate models, Model-A and Model-B, with corresponding geometries of the caustics shown in Figure 3 and 4. The non-parallax, parallax and parallax plus lens orbital motion (LOM) models are plotted in Figure 2 and denoted by the unmarked model names, the extra character P and P+LOM, respectively.Subscripts "+" and "−" refer to degenerate parallax models with corresponding signs of the impact parameter u 0 .Parameters of the models are listed in Tables 2 and 3.
The event timescales of both Model-A and Model-B, t E ∼ 70 days and 110 days respectively, are longer than the typical timescale of t E ∼ 30 days for events toward the Galactic bulge (Mróz et al. 2019), from which we anticipate a heavy lens.The Model-A resembles the Model-B in the binary lens parameters q, s and α, but differs in t E and ρ.The larger t E value and the smaller ρ value of Model-B implies a larger θ E compared to that of Model-A.
The non-parallax model of Model-A already fits well to the light curve whereas Model-B shows a deviation around HJD ′ ∼ 6760 due to cusp re-approach after the caustic exit, resulting in ∆χ 2 ∼ 45 between Model-A and Model-B.Once the higher order effects are considered, Model-B significantly improves from the nonparallax model by avoiding the cusp re-approach and ∆χ 2 ∼ −60 and ∼ −85 for the parallax and parallax plus LOM models.In contrast, Model-A only improves by ∆χ 2 = −6, showing no distinct evidence of the parallax signal.Therefore, we conclude a higher order effect of either parallax, lens orbital motion or xallarap needs to be involved in Model-B but not in Model-A.
In order to take the higher order effect into account for Model-B, we begin with the parallax alone model that already provides a good fit to the flat feature around HJD ′ ∼ 6760 and then add LOM for getting more rigorous parallax parameter values and errors.Our modeling suffers a degeneracy known as the ecliptic degeneracy, where light curves with the parameters (u 0 , α, π E,N , γ ⊥ ) = −(u 0 , α, π E,N , γ ⊥ ) appear nearly identical if direction of the Sun's acceleration is constant (Skowron et al. 2011).This occurs when a deviation due to parallax is short timescale or the event coordinate is close to the ecliptic plane which is true for the Galactic bulge.We identify the degenerate models with the subscripts "+" and "−" taken from the sign of u 0 .The "+" and "−" models have similar parameter values except the sign and are indeed indistinguishable from the light curve and χ 2 .The inclusion of LOM improves the fitting by ∆χ 2 ∼ −25 and reduces uncertainty of the parallax vector from ∼ 30% to ∼ 15% in the positive π E,N regime and ∼ 6% to ∼ 2% in the negative π E,N regime, as can be noticed in Tables 2 and 3 and Figure 5 where we show the posterior distribution of the parallax vector components with uniform prior of π E,N and π E,E .However, we find the ∆χ 2 mostly comes from the baseline systematics but not from the magnification part of the light curve, which indicates that an effective χ 2 improvement due to the microlensing signal is negligibly small.We also examined the xallarap effect by simultaneously fitting with the parallax parameters, but it changes neither the model nor the parallax vector much.We note that values of t E and ρ are consistent between the parallax alone model and the parallax plus LOM/xallrap model.
Despite the strict constraint of the parallax parameters imposed by simultaneous fitting of the parallax ef- fect and LOM, we confirmed that the deviation around HJD ′ ∼ 6760 is also explainable by the LOM or xallarap alone, which accepts null detection of the parallax.We thus conclude the parallax parameters cannot be determined contrary to the small model uncertainty.
There is also a large dispersion in the distribution of ρ as shown in Figure 6 where the posterior distributions of ρ are plotted.The dispersion is due to the sparsity of data points taken during the caustic entry and exit.The distributions are well characterized by the superposition of two Gaussian distributions, and the best-fit finite source parameter of the models lies either of the bimodal peaks.This bimodal feature is associated with the uncertainty of ρ, causing it to be larger.We calibrated the source magnitude in the instrumental OGLE I -and V -band to the standard Kron-Cousins and Johnson system using the following equations from Udalski et al. (2015): where i DB and v DB are the instrumental I and V , µ = (1 − ϵ V + ϵ I ) −1 , ∆ZP I = 0.018, ∆ZP V = 0.158, ϵ I = −0.005± 0.003 and ϵ V = −0.077± 0.001 for the field of OGLE-2014-BLG-0221.The calibrated color and magnitude (V − I, I) S are plotted onto the colormagnitude diagram (CMD) of the stars within 2 ′ of the event coordinate (Figure 7).The central color and magnitude of the red clump giants (RCGs) population marked in the CMD is representative of the bulge RCGs.We followed the standard procedure adopted in Yoo et al. (2004) to correct the effect of reddening and extinction due to interstellar dust, under the assumption that the source experiences the same reddening and extinction as the bulge RCGs.The centroid of the RCGs was (V − I, I) RCG = (2.358,15.947)±(0.009,0.023) from the color-magnitude distribution, and the extinction-free color and magnitude of the bulge RCGs were known as (V −I, I) RCG,0 = (1.060,14.349) ± (0.070, 0.040) toward the event coordinate (Bensby et al. 2011;Nataf et al. 2013), and for which the reddening and extinction were estimated as (E(V − I), A(I)) = (1.298,1.598) ± (0.071, 0.046).Applying the same reddening and extinction and using the earlier empirical relation, we found the intrinsic color and magnitude of the source (V − I, I) S,0 and the angular source size θ * for Model-A, Model-B P± and Model-B P+LOM±, that let us calculate θ E and µ rel .Those values are given in Table 4.

LENS PROPERTY
Since parallax could not be measured accurately, we used a Bayesian approach to infer the lens physical properties.The microlensing event simulation code (Koshimoto & Ranc 2021) with a parametric Galactic model toward the Galactic bulge developed by Koshimoto et al. (2021) was used to generate artificial microlensing events and obtain posterior distributions of the physical parameters.Likelihood of the measured parameters, t E and θ E , was included as the observed constraints.An upper limit of the lens brightness was also set as I lim = 18 mag, determined based on the blending magnitude but about 2 mag brighter to make the analysis conservative.From the OGLE-III catalog, we also confirmed there is no potential source within 2 ′′ of the event coordinate that is brighter than the limit.Apparent lens magnitudes of the simulated events were estimated using the massluminosity and color-color relations of main sequence stars (Kroupa et al. 1993;Kenyon & Hartmann 1995) and the extinction law (Nishiyama et al. 2009).The following extinction model was assumed from Bennett et al. (2015) for dependence on the lens distance, where A L is the extinction experienced by the lens and A RC is the extinction of the bulge RCGs toward the event coordinate.The equation presumes dust in foreground of the lens with a dust scale height of 0.1 kpc.
In addition to a stellar luminous lens, we also generated a remnant dark lens in accordance with the PARSEC isochrone models (Bressan et al. 2012;Chen et al. 2014;Tang et al. 2014) of stellar evolution and the remnant initial-final mass relation (Lam et al. 2020) implemented in the simulation code.Natal kick velocities of 350 and 100 km/s are assumed for neutron stars and black holes respectively, following Lam et al. (2020).Figures 8, 9 and 10 show the posterior probability density distributions of the physical parameters for Model-A, Model-B P + and Model-B P+LOM + with 10 6 simulated microlensing samples accepted under the constraints, and Table 5 lists the physical parameters.We did not simulate events for Model-B P − and Model-B P+LOM − since their likelihood and constraints on the Bayesian analysis are practically same as that for Model-B P + and Model-B P+LOM + , respectively.
The event simulation results assuming Model-A indicate the lens system is most likely a planetary system consisting of a late-type star orbited by a gas giant.The simulation also generated remnant samples; white dwarf lenses at almost identical parameter ranges as main sequence lenses occupy ∼ 70% of the remnant distributions.A bimodal feature appears in the lens distance distribution because the prior probability is weighted more to the Galactic center region whereas the model requires a nearby lens for a low mass system, i.e. main sequence and white dwarf lenses, to explain the longer t E and larger θ E , which are proportional to the square root of π rel = (1 au)(1/D L − 1/D S ).On the other hand for Model-B, the result supports a conclusion that the lens system is most likely composed of a remnant orbited by a gas giant, brown dwarf or red dwarf depending on the host mass, which is a reasonable consequence of the longer t E and larger θ E likelihood compared to that suggested by Model-A.Especially, it is notable that the predicted percentage of remnant lenses is 100% for Model-B P+LOM + with more than 95% of the distributions occupied by black hole lenses.Discrepancies between the distributions of the Model-B P + and Model-B P+LOM + parameters as well as fractions in the lens types mostly comes from the broad likelihood distribution of the Model-B P + parameters due to the dispersion in the finite source effect parameter.Upper end of the main sequence distributions are constrained by the upper limit imposed on the lens brightness.

SUMMARY AND DISCUSSION
The results of the light curve modeling and investigation of the source and lens properties were shown in the preceding sections.We found two degenerate models, Model-A and Model-B, that have similar companion parameter values of q ∼ 5×10 −3 and s ∼ 1.1 but dissimilar microlensing parameter values of t E ∼ (70, 110) days and ρ ∼ (5, 1) × 10 −4 .Due to the long event timescale and large angular Einstein radius, both models favored a nearby heavy lens solution; furthermore, a Bayesian analysis including remnant populations in the lens system prior revealed the lens to be a remnant candidate.

Interpretation of the Lens System
Although several combinations of the host and companion object types are proposed regarding the posterior distributions as shown in Table 5, they can be divided roughly into three pairs, a gas giant planet with a main sequence star, a gas giant with a remnant, and a brown dwarf or red dwarf with a black hole (BH).
Both close-in and distant giant planets around main sequence stars are found to date by various survey techniques and statistically studied for their occurrence around different spectral types of hosts.For microlensing, Suzuki et al. (2018) compared the statistical analysis result of planet occurrence presented in Suzuki et al. (2016) to the core accretion model (Ida & Lin 2004;Mordasini et al. 2009) and confirmed excess in microlensing planets beyond q ∼ 10 −4 , including a factor ∼ 5 discrepancy for 10 −3 < q ≤ 0.03.The favored value for q for event OGLE-2014-BLG-0221 will add to that discrepancy.The MOA collaboration is preparing an extended analysis of Suzuki et al. (2016) with extended data beyond 2007 − 2012 in which OGLE-2014-BLG-0221 will be included.There is other observational evidence from various observational methods that support some modification to the traditional planet formation scenarios of giant planets.A recent study of transiting planets discovered by TESS indicates a higher occurrence rate of giant planets around low mass stars than the rate predicted by the core accretion model (Bryant et al. 2023).The discovery of four distant giant planets around HR 8799 by direct imaging suggests that gravitational instability (Boss 1997) would play an important role in a planet formation (Marois et al. 2008(Marois et al. , 2010)).
There are a few cases of exoplanets discovered around a white dwarf (WD) or a neutron star, including the first exoplanets detected in 1992 around a pulsar, PSR B1257+12 (Wolszczan & Frail 1992).However, most of them represent extreme environment; for example, PSR B1620-26AB b is a circumbinary planet around a pulsar-WD binary (Thorsett et al. 1993;Sigurdsson et al. 2003), WD 0806-661 b has a very large separation of 2500 au from its host (Luhman et al. 2011) and WD1856+534 b has a small separation of 0.02 au ( Vanderburg et al. 2020).Only MOA-2010-BLG-477Lb is a currently known planet orbiting about a few au away from a WD (Blackman et al. 2021).It is still uncertain how a solar analogous planetary system evolves along with its host star.Several mechanisms are proposed such as common envelop during the giant phase (Paczynski 1976) that results in short planetary orbits or stellar mass loss that pushes planets outward (Veras 2016).Second generation exoplanets, those formed during the post main sequence phase, are also proposed as an alternative scenario for some of the observed systems (Perets 2010;Ledda et al. 2023).Despite its importance of finding more observational samples, observational bias hinders the solid detection of a WD planetary system.The low luminosity of WDs makes astrometry and transit monitoring difficult, the non-characteristic features of a WD spectrum prevents measuring radial velocity, and direct imaging is biased toward wide orbits as WD 0806-661 b shows.In contrast, microlensing has an advantage as it does not rely on host star brightness, and a microlensing survey toward the Galactic center has a peak sensitivity of planet-host separation at a few au, that is suitable for filling the gap between close-in to distant planets around WDs. Microlensing can also detect a planet around a X-ray quiet neutron star with its capability of finding dark objects that is not achievable by any other method.Although the pulsar timing method is the only method successful to date for detect-ing a planet around a neutron star 2 , microlensing would shed light on a cold planet around an unseen neutron star if the host of OGLE-2014-BLG-0221 is a neutron star.
More than a dozen stellar mass BH binaries are known but are mostly found by light from X-ray binaries where accretion of material from a companion occurs and is therefore a closely packed system (e.g.Corral-Santana et al. 2016).A few non-interacting black hole candidates are also identified by radial velocity (e.g.Shenar et al. 2022;Mahy et al. 2022) and astrometry (El-Badry et al. 2023), in which ∼ 1.4 au is the widest orbital separation between the candidate and its companion.Several population synthesis studies predict that large fractions of BH luminous companion binaries should have a wide separation such that orbital period becomes more than years (e.g.Chawla et al. 2022), much longer than any of the confirmed candidates.We estimated OGLE-2014-BLG-0221 has a projected separation of ∼ 10 au for the case of a black hole host; thus, it would represent the longest BH binary separation ever known and belong to the theoretically predicted population of wide orbit BH binaries.Moreover, the companion is supposed to be a very low mass star or a brown dwarf, and in either case is too faint to be detected by photometry, indicating such a system is unlikely found by any other method besides microlensing.

Future Follow-up Observation
Once the source and lens are separated enough following their proper motions after several years of the event peak, we would be able to resolve the source and lens with high resolution imaging (e.g.Hubble Space Telescope, Keck adaptive optics).This idea was first developed in Bennett et al. (2006Bennett et al. ( , 2007)), and the number of identifications has increased in recent years (e.g.Bennett et al. 2015;Bhattacharya et al. 2018;Terry et al. 2021).Measurement of the lens-source separation and the lens brightness allows strong constraints to be placed on the mass-distance relation of the lens.Here, we consider the detectability of the lens using the future high resolution imaging of OGLE-2014-BLG-0221.
Figure 11 is the posterior distribution of the predicted source-lens separation for each degenerate model with the parallax parameters included, derived from the MCMC parameter chains.The separation is computed from the heliocentric proper motion vector: 2 7 planets are confirmed including 3 plants belong to the PSR B1257+12 system (Akeson et al. 2013).where πE is the unit vector of the microlensing parallax and v ⊕ is the Earth's projected velocity at t 0 .While the direction of the separation is largely dependent on the microlensing parallax vector, the conversion from the geocentric to the heliocentric reference frame does not change the separation itself much unless the relative parallax is too large.We estimated the current separation is 41 ± 7 and 45 ± 14 mas for Model-A P ± , 103 ± 28 and 98 ± 13 mas for Model-B P ± , and 139 ± 20 and 132 ± 28 mas for Model-B P+LOM ± , respectively.Hence, we expect the source and lens are separated enough to conduct the high resolution imaging, otherwise the separation is unexpectedly small owing to its uncertainty.By measuring the current position, we would discriminate the models as the direction and magnitude of the separation are likely different among the models.The large expected source-lens separation of Model-B would be due to the small lens distance or the kick velocity caused when a NS or BH forms as discussed in Lam et al. (2020).
The relative source-lens brightness is also a major concern when resolving the source and lens.Similar brightness is preferred to identify both source and lens; however, Bhattacharya et al. (2021) demonstrated that the lens can be identified even if it is a few magnitudes fainter than the source.For the case of OGLE-2014-BLG-0221, photometric detection of the lens indicates that the system most likely follows Model-A.From the Bayesian posterior distribution, we obtained the es-timation of the apparent lens magnitude (H L , K L ) = (17.9+1.1 −1.1 , 17.7 +1.0 −1.1 ) for Model-A.This should be luminous enough as, for example, Blackman et al. (2021) determined the detection limit of their high resolution Keck image as 21.1 mag in H -band.In comparison to the lens brightness, the similar apparent source magnitude is expected from the modeling as (H S , K S ) = (18.1,17.9) ± (0.8, 0.8) for Model-A, making the event ideal for the high resolution imaging.Furthermore, Model-B measures the different apparent source magnitude, (H S , K S ) = (18.9,18.7) ± (0.8, 0.8) and (18.8, 18.6) ± (1.2, 1.2) for Model-B P + and Model-B P+LOM + , from that of Model-A.This indicates the model degeneracy would be disentangled with observations of the source even though the lens is not observable for the case of a remnant lens.Future high resolution imaging is highly important for characterizing OGLE-2014-BLG-0221.

Figure 1 .
Figure 1.The light curve of OGLE-2014-BLG-0221.The black, red and green data points are the reduced photometric data from the MOA-Red, OGLE I -and V -band observations, respectively.The blue and green lines represent the best-fit binary lens single source model and best-fit single lens single source model.The lower panel shows the residuals from the best-fit binary lens model.

Figure 2 .
Figure 2. Model light curves of the two degenerate models around the anomaly.The Model-A, Model-A P±, Model-B, Model-B P± and Model-B P+LOM± are plotted.Since different source and blending fluxes between the Model-A and Model-B make different magnification, magnification of the Model-B is scaled by source and blending fluxes of the Model-A for comparison.The lower residual panel shows the data points binned in 1 day interval.

Figure 3 .
Figure 3. Caustic geometry corresponding to Model-A for versions with non-parallax and parallax (P±).The caustic is shown in the red curved line.The blue line is the source trajectory.

Figure 4 .
Figure 4. Caustic geometry corresponding to Model-B for versions with non-parallax, parallax (P±) and parallax plus lens orbital motion (P+LOM±).The caustics at HJD ′ = 6728.6,6738, 6746 and 6760 are plotted for P+LOM± with the instantaneous source positions in the red open circles whose sizes are not in scale.

Figure 5 .
Figure 5. Posterior distribution sampled by the MCMC sampler.Colors represent ∆χ 2 from the lowest χ 2 parallax plus LOM model.

Figure 6 .
Figure 6.Probability density functions (PDFs) of the normalized angular source radius ρ from our MCMC chains for (a) Model-A, (b) Model-B P+ and (c) Model-B P+LOM+.The bimodal distribution of each model except Model-B P+LOM+ is approximated by the Gaussian Mixture Model, and two components of the best fit is plotted in orange and green lines.
4. SOURCE PROPERTYOnce the finite source effect is detected, the angular Einstein radius θ E = θ * ρ can be determined by measuring the angular size of the source, and the relative proper motion µ rel = θE tE can also be estimated, combining with t E measured from the light curve.An empirical relation of the angular source size as a function of the cousins I and V band is derived based on the result ofBoyajian et al. (2014).By restricting to stars with 3900 K < T eff < 7000 K, optimal to FGK stars,Fukui et al. (2015) finds: log[2θ * /mas] = 0.5014 + 0.4197(V − I) − 0.2I.(5)

Figure 7 .
Figure 7. (V − I, I) Color-Magnitude Diagram of stars in the OGLE-III catalog (Szymański et al. 2011) and the Hubble Space Telescope catalog (Holtzman et al. 1998).The black dots are the OGLE-III catalog stars located within 2 ′ of the event OGLE-2014-BLG-0221, and the green dots are the Hubble Space Telescope catalog stars in a field of Baade's Window, whose colors and magnitudes are shifted to match the red clump giant centroid of the OGLE-III catalog.The red spot indicates position of the red clump giant centroid, and the orange, blue, purple, cyan and magenta spots are source positions of Model-A, Model-B P± and Model-B P+LOM±, respectively.

Figure 8 .
Figure 8. Posterior probability density distributions of Model-A.The distributions are classified into main sequence stars, white dwarfs, neutrons stars and black holes and separated by colors accordingly.In the last plot, lens brightness of the main sequence samples is shown where the vertical blue line is the median value, and the cyan and purple regions represent 68.3 and 95.4% credible interval of the distribution.The apparent source brightness is also plotted by the median magenta line with the magenta dashed line indicating the 68.3% interval.

Figure 10 .
Figure 10.Same as Figure 8 for Model-B P+LOM+.The plot of the lens brightness is not shown as no main sequence sample has been generated in the simulation.

Figure 11 .
Figure 11.Current source-lens relative position.The distributions of Model-A P±, Model-B P± and Model-B P+LOM± are plotted.The contour from dark to light shows 39.3, 86.5 and 98.9% highest density regions.

Table 1 .
Data Set and Rescaling Parameters

Table 2 .
Model-A Parameters Note-Best-fit values with uncertainties corresponding to the 68% credible intervals of the MCMC posterior distributions.

Table 3 .
Model-B Parameters

Table 5 .
Estimated Lens Physical Parameters from the Bayesian Analysis