© The Authors 2025

1 Introduction

More than three decades have passed since Mao & Paczynski (1991) first suggested that extrasolar planets could be detected in microlensing events. During this time, over 230 planets¹ have been discovered using this method. This number may seem small compared to discoveries made by other major planet-detection techniques, but the strength of microlensing lies in its unique sensitivity to low-mass and distant planets. Instead of analyzing the light of the planet’s host star, microlensing uses the light of an unrelated background star to probe the foreground planetary system. This makes the technique sensitive to, for example, planets at and beyond the snow line of their host stars, as well as to those located at distances from the Galactic center that are currently inaccessible by any other method (Tsapras 2018).

Another advantage of microlensing is its usefulness in demographic studies². In microlensing, planetary occurrence rates are typically derived as a function of mass ratios (q) and projected separations in Einstein ring radius units (s), calculated using a set of events and the detection efficiency of the survey that detected them (Mróz & Poleski 2024). However, these analyses are affected by uncertainties and biases that require significant effort to disentangle (e.g., Gould et al. 2010; Udalski et al. 2018). These factors are inherent to the sample of analyzed microlensing events, starting with observational data collected by numerous surveys using various instruments and under differing conditions. Furthermore, variations in data reduction and selection processes, along with subjective human factors such as publication bias (Yang et al. 2020), contribute to these uncertainties. Collecting a sample of events detected from homogeneous data gathered by a single instrument greatly minimizes these issues.

We present the discovery of the planetary microlensing event OGLE-2015-BLG-1609. The planetary signal in this event was identified solely through a single survey, making it well-suited for the statistical analysis of microlensing planets. The analysis turned out to be more complex than typical, due to low-level systematic trends in the photometry. We find that there are three different topologies that could explain the light curve of the event; two of them have very similar Bayesian evidence, and hence, we were not able to distinguish between them. We estimated the mass of the planet host, which overlaps with masses of both stars and brown dwarfs.

One of the most detailed studies on the planetary occurrence rate in microlensing events was done by the Microlensing Observations in Astrophysics (MOA) collaboration, using their 2007–2012 sample (Suzuki et al. 2016). This sample comprised 1474 events, including 23 planetary events. Suzuki et al. (2016) find that the planetary occurrence rate follows a broken power law of q. However, due to the lack of events in the sample with q < 10^−4.5, which is close to the break in the power law (q_br), fitting q_br led to a high uncertainty for all the parameters. Therefore, the authors fixed this value to q_br = 1.7 × 10⁻⁴ and concluded that beyond the snow line, most planets should have Neptune-like masses. Later, known planets with q < 10⁻⁴ were analyzed by Udalski et al. (2018). Combining planets detected using both survey and follow-up data forced authors to apply an alternative inference method (the V/V_max method; Schmidt 1968). Udalski et al. (2018) confirmed the break in the mass ratio power law but at higher masses, q_br = 2 × 10⁻⁴. The break was further analyzed by Jung et al. (2019) using a sample of 15 planets with q < 3 × 10⁻⁴ and assuming that planet-detection sensitivity as a function of q can be approximated by a simple power law. Their results suggest that the break is at a much smaller mass ratio (q_br ≈ 0.55 × 10⁻⁴) and the slope of the distribution at lower mass ratios is much steeper compared to the findings by Suzuki et al. (2016).

Significant efforts are currently being undertaken, using the Korea Microlensing Telescope Network’s (KMTNet) semiautomated algorithm Anomalyfinder (Zang et al. 2021), to build a large, uniformly selected sample of planetary events. So far, the algorithm has identified about 100 planets in KMTNet’s observations from the years 2016–2019 (Zang et al. 2021; Hwang et al. 2022; Wang et al. 2022; Gould et al. 2022; Zang et al. 2022; Jung et al. 2022; Zang et al. 2023; Jung et al. 2023; Shin et al. 2023; Ryu et al. 2024; Shin et al. 2024), including a few with q < 0.5 × 10⁻⁴. A future statistical analysis of this sample will significantly improve constraints on planetary occurrence rates, in particular in the low-mass regime. The projected separation and the mass ratio are not the only parameters that can be derived from microlensing events. By measuring the so-called secondary effects in the light curve, it is possible to obtain absolute masses, the distance to the system, and the projected separation in physical units. For this purpose, a pair of finite-source and microlensing parallax effects is typically used. Despite the fact that both of these effects are measured only in ~20% of planetary events, some initiatives were taken to measure how the distribution of the planets depends on the host mass and Galactocentric distance (M_h, R). The hypothesis of an equal abundance of planets in the bulge and the disk was tested by Penny et al. (2016) using a sample of 21 systems. However, the small sample of the parallax measurements, hindered by systematic uncertainties, prevented them from reaching a definitive conclusion. Additionally, some nearby planets (< 2 kpc) in the sample have now been shown, through analysis of high-angular-resolution images, to be at larger distances than originally reported, for example, MOA-2007-BLG-192 (Terry et al. 2024). Using a slightly larger sample (28 planets), Koshimoto et al. (2021b) derived the planetary occurrence rate as a function of both host mass (M_h) and Galactocentric distance (R), $P_{\text{host\hspace{0.17em}}}\propto M_{\mathrm{h}}^m{R}^r$, by comparing the observed distribution of lens-source relative proper motion (μ_rel) to Galactic model predictions for a given Einstein radius crossing time (t_E). Their results suggest that the probability of hosting a planet increases with the Galactocentric distance, but this dependence is not statistically significant, with r = 0.2 ± 0.4. The uncertainties in the exponents of M_h and R were reduced by Nunota et al. (2024) by utilizing model distributions of both μ_rel and t_E and dividing systems into two subsamples, with mass ratios below and above q = 10⁻³. Their analyses indicate that beyond the snow line, massive planets are more likely to be around more massive stars, while the frequency of low-mass planets does not depend on the mass of the host.

The planetary event analyzed here will be included in the statistical studies of MOA or Optical Gravitational Lensing Experiment (OGLE) survey planets. Such studies can shed more light on the detailed aspects of the planetary occurrence rate.

This paper is organized as follows: in Sect. 2 we describe the observational data and their preliminary analysis. In Sect. 3, we detail the model fitting to the OGLE observations alone. In Sect. 4, we present the final results of modeling using all available datasets and the process of estimating the physical parameters of the system. Finally, we summarize our results in Sect. 5.

2 Observational data

2.1 Collection

The microlensing event OGLE-2015-BLG-1609 was first detected by OGLE on 2015 June 16, heliocentric Julian date (HJD) ~2457219. The detection was made by the Early Warning System (EWS³; Udalski 2003) through the analysis of data collected by the 1.3-m Warsaw Telescope, equipped with the 1.4 deg² field of view mosaic CCD camera at the Las Campanas Observatory in Chile (Udalski et al. 2015). The event was located at (RA, Dec)₂₀₀₀ = (18:03:17.71, -26:54:25.6) in equatorial coordinates, or (l, b) = ($3_{\cdot }^{\circ }72,-2_{\cdot }^{\circ }34$) in Galactic coordinates, which is inside the BLG511 field of the OGLE-IV survey. This field was one of the OGLE high-cadence Galactic Bulge fields, observed four times per night on average. OGLE images were primarily taken in the I band, with occasional observations in the V band conducted approximately once every few days. The reduction of the OGLE data was done using a variant of difference image analysis (DIA; Crotts & Tomaney 1996; Alard & Lupton 1998) optimized by Wozniak (2000) and Udalski et al. (2008).

The event was also observed by MOA’s 1.8 m telescope at Mt. John University Observatory in New Zealand. It was independently alerted on 2015 June 30, HJD ~ 2457233, as MOA-2015-BLG-412 on the MOA alerts web page⁴. MOA observations were primarily done in the custom wide MOA-Red filter, approximately equal to the sum of R− and I− band filters. MOA data reductions were performed using a variation of DIA presented by Bond et al. (2001). Additionally, MOA observed in the V band, but these few observations were excluded in this analysis due to large error bars.

On 2015 September 3, HJD ~ 2457269, the light curve visibly deviated from the standard Paczynski (1986) model, indicating the presence of a planetary anomaly. This triggered follow-up observations on three 1 m Las Cumbres Observatory (LCO) telescopes located at the South African Astronomical Observatory, operated by the RoboNet collaboration (Tsapras et al. 2009), and on the 1.54 m Danish Telescope at La Silla Observatory, Chile, operated by the Microlensing Network for the Detection of Small Terrestrial Exoplanets⁵ (MiNDSTEp) collaboration (Dominik et al. 2010). RoboNet data reductions were performed using a modified DIA algorithm (DANDIA; Bramich 2008). Observations on the Danish Telescope were collected using a high frame-rate electron-multiplying CCD (Skottfelt et al. 2015) and were reduced using a variation of the DANDIA algorithm optimized for this instrument. The light curve of the event is plotted in Fig. 1, with a closer view of the planetary anomaly provided in Fig. 2. Some data were collected using other telescopes but due to poor coverage of the event, large error bars, or the lack of reliable reductions, we excluded them from this analysis.

Fig. 1 Light curves of the microlensing event OGLE-2015-BLG-1609 with three microlensing model topologies for positive values of the impact parameter, u0. The difference between the topologies is visible only during the planetary anomaly (HJD ~ 2457270; see Fig. 2). The photometry comes from the OGLE, MOA, RoboNet, and MiNDSTEp projects.

Fig. 2 Same as Fig. 1 but for the light curve of the planetary anomaly in OGLE-2015-BLG-1609.

2.2 Preliminary analysis

In the initial analysis of OGLE-2015-BLG-1609 based solely on the OGLE survey data, we faced several challenges. In models that include the microlensing parallax effect, we obtained well-constrained but astrophysically unlikely values for the northern component of the parallax vector (π_E,N), which should be close to zero according to the Galactic models (e.g., Lam et al. 2020; Koshimoto et al. 2021a). Additionally, at the end of the event and before the annual break in the visibility of the Galactic bulge, a trend at HJD 2457300–2457335 was noticed in the residuals from the best fit. We tried to solve these issues by including additional effects in microlensing models. First considered was the binary source effect (2L2S). In this model, the brightening at the end of the event was caused by the presence of a second star in the source system. Inclusion of this effect improved χ² of the fit, but did not change the value of π_E,N. Since only half of the second source brightening was covered by data points, we conclude that this trend is more likely an end-of-the-season effect than an astrophysical one, and we rule out the binary source explanation.

We also considered additional changes in the geometry of the event caused by either the orbital motion of the lenses or the orbital motion of the source (the xallarap effect). In the case of lens orbital motion, we did not achieve reasonable constraints on the effect’s parameters, nor did we improve the χ² statistics of the fits. This could be attributed to the lack of distinctive features in the light curve, such as caustic structures, that can serve as time references (An & Gould 2001). Similarly, in models with the xallarap effect, there was no improvement in the χ² statistics. Additionally, the resulting values of the orbital period of the source aligned with spurious periods of low-level amplitude recurring in the OGLE photometric data, as found in previous studies (Mróz et al. 2023). These periods are associated with the structure of the data itself. For this reason, we do not consider the results of xallarap models to be trustworthy. Including lens orbital motion or xallarap effects in the models did not change the values of π_E,N. The overestimation of the parallax vector may be caused by other unassociated photometric variability. To address the potential trends in our data we incorporated Gaussian processes into modeling of the event. This approach has been successful in modeling the variability of the microlensing source, demonstrated in the analysis of the event OGLE-2017-BLG-1186 (Li et al. 2019). Unfortunately, it led to an increase in the π_E,N value with even smaller uncertainties. Simultaneous modeling both the variability of the source star and the parallax effect requires careful selection of the baseline time range. A long baseline is needed to reliably constrain the long-term stellar variability. However, extending the baseline much beyond the duration of the microlensing event can adversely affect parallax vector measurements by introducing uncorrelated trends.

Finally, to ensure the absence of long-term observational trends in the light curve of the event, we analyzed observations of the nearby stars in the OGLE database. We selected 16 constant stars in the OGLE database around 30″ from the source star, with similar brightness and color. We arranged their observations in week-long bins, then calculated the average between all the chosen stars. The resulting trend did not exceed the average magnitude uncertainty of the analyzed data and thus could not affect the microlensing modeling.

After trying all the approaches described above, we decided to include the data collected by the MOA group in our analysis. Models fitted to the MOA data alone and to the combined MOA and OGLE datasets provided more typical results for the π_E,N value, which was consistent with zero within the 2σ limit. This led us to make a step back and take a closer look at the reductions of the OGLE images. In the DIA reduction method, photometric measurements are obtained by subtracting a reference image from the acquired images. The reference images are created by averaging a selected dozen or so images taken under the best seeing conditions. For the initial OGLE reductions, reference images were selected from those taken at the start of the 2010 season. Three years after the event, in 2018, the OGLE team prepared new reference images. Since they were closer in time to the analyzed event, we used them to re-extract the photometry. We fitted the model with parallax effect to this new OGLE dataset. This time we obtained a distribution of π_E,N that was consistent with zero within the 2σ limit and, therefore, in a good agreement with the Galactic model. The photometry extracted using DIA can be affected at a low level by the proper motion of the source star. Hence, changes in the position of the source star can influence microlensing parallax measurements. We present in Table 1 proper motion measurements from the Gaia DR3 (Gaia Collaboration 2023) and the internal OGLE time-series astrometry, which has been matched to the Gaia DR3 reference frame. We note that the uncertainties in the OGLE measurements are purely statistical and likely underestimated. The value of the renormalized unit weight error (RUWE > 1.5; Table 1) for the Gaia DR3 asymmetric solution indicates that the parallax measurements from Gaia DR3 cannot be attributed to either the lens or the source in the event (Wyrzykowski et al. 2023).

In all cases, the binary-lens models fitted to individual datasets from both the “old” and “new” OGLE reductions, as well as MOA observations, have significantly lower-χ² statistics compared to the single-lens models fitted to the same datasets. This enables the inclusion of this event in the planetary occurrence rate studies of both surveys. The comparison is shown in the Table 2. In the subsequent investigation, we used new OGLE reductions only.

2.3 Data preparation and photometric uncertainties

We performed the color calibration of the OGLE I – and V-and measurements, so that magnitudes reported in this work are in the standard I (Cousins) and V (Johnson) pass-bands. To clean datasets, we excluded all observations with error bars larger than five times the median of the error bars of nearby points. Additionally, using one of the early models of the event, we removed 3σ outliers from the fitted light curve.

Photometric pipelines sometimes struggle with estimating the impact of systematics in the data; uncertainties in the raw measurements are often underestimated and have broader tails than a Gaussian. To renormalize photometric error bars, we scaled them using the formula for the i-th data point (Yee et al. 2012): $${\sigma }_{\mathrm{n}\mathrm{e}\mathrm{w},i}=\sqrt{{\left(k{\sigma }_i\right)}^2+e_{min}^2},$$(1)

where σ_i is an initial error bar, k and e_min are re-normalization coefficients. Since we used new reference images for OGLE data reductions, the values of k and e_min from the empirical model for this dataset reported by Skowron et al. (2016) could be inaccurate. For the sake of coherent treatment of each dataset, we regarded k and e_min for each dataset as additional parameters of the model and sampled with Markov chain Monte Carlo (MCMC). We assumed that error bars should be Gaussian in flux space, and used the modified likelihood function in the form of $$\begin{array}{rl}\ln \mathcal{L}& =-\frac{1}{2}\sum _{i=1}^N[{\left(\frac{f_i-f_{\text{mod},i}}{{\sigma }_{\text{new},i}}\right)}^2+\ln \left(2\pi {\sigma }_{\text{new},i}^2\right)]\\ & =-\frac{1}{2}[{\chi }^2+\sum _{i=1}^N\ln \left(2\pi {\sigma }_{\text{new},i}^2\right)],\end{array}$$(2)

where f_i − f_mod,i represents the difference between the measured and modeled flux values at the i-th observation. All the modeling presented in this work was initially done using this approach. Then, the error bars scaling parameters were fixed to the values from the models with the lowest χ², and the modeling was repeated.

3 Single instrument data modeling: OGLE

As the planetary anomaly is well sampled in the OGLE observations alone, initially we focused on just this dataset in our analysis. The observed brightening in the OGLE-2015-BLG1609 light curve, except for the anomaly, can be described well by the classical Paczynski (1986) model. In this single-lens point-source model (1L1S), microlensing magnification is described by $$A(u)=\frac{u^2+2}{u\sqrt{u^2+4}},$$(3)

where u(t) denotes the projected separation between the source and the lens at a given time, scaled to the angular Einstein radius (θ_E). This separation is related to the impact parameter of the lens-source approach (u₀) and the time of closest approach (t₀) through the relation $$u(t)=\sqrt{u_0^2+\tau (t{)}^2}\text{\hspace{0.17em}and\hspace{0.17em}}\tau (t)=\frac{t-t_0}{t_{\mathrm{E}}}.$$(4)

Taking into account that not all the measured light is magnified, the flux change as a function of time is given by $$f(t)=f_{\mathrm{S}}A+f_{\mathrm{B}},$$(5)

where f_S and f_B are the fluxes of the source and blended objects, respectively. As the starting point of our analysis, we excluded data points collected during the anomaly (HJD 2457267–2457273) and fitted the Paczyński curve to the remaining observations. The 1L1S parameters (t_E, t₀, and u₀) were used as initial values for further modeling. Next, we considered the binarylens single-source (2L1S) model, which is outlined in Sect. 3.1. We extended this model by including the microlensing parallax effect (Sect. 3.2). Moreover, due to unreliable parameter constraints, we also imposed priors on π_E based on a Galactic model, which is described in Sect. 3.3.

All the microlensing modeling in this work was done using the MulensModel Python-code (Poleski & Yee 2019), with the affine invariant MCMC ensemble sampler (Goodman & Weare 2010) implemented in the emcee code (Foreman-Mackey et al. 2013). The number of chains and steps for each MCMC modeling run was adjusted to achieve the convergence. All runs with fixed error bars scaling parameters had a mean acceptance fraction between 0.287 and 0.458, the number of steps is at least 34 times longer than the mean autocorrelation time, and no thinning was applied (for definitions, see Foreman-Mackey et al. 2013).

3.1 Binary-lens model

The planetary anomaly in the light curve led us to consider a binary-lens model with finite-source effect (2L1S). This model requires extending the number of parameters from three in Paczyński’s model to seven. Those additional parameters are typically: s – the projected separation between the primary lens and the companion, q – the ratio of the companion’s mass to the primary lens mass, α – the angle between the source trajectory and the binary axis, and the parameter of the finite source effect: ρ – the radius of the source in the Einstein radius units. In our modeling, we used this parametrization. Other alternative parametrization were not suitable, as most of them are connected to features of the light curve involving caustic crossings, which are absent in the analyzed event.

To find all the degenerate solutions, we conducted a grid search in the log s − log q parameters space. We used fixed values of s evenly spaced between −0.15 < log s < 0.35 and q spaced between −4.5 < log q < −2, and optimized with MCMC method values of other parameters (t_E, t₀ u₀, ρ, and α). As the absence of a visible caustic crossing results in poor constraints on α, we split our modeling into two steps. First, we spread uniformly starting values of α and ran EMCEE with 300 walkers. After confirming that all models across log s − log q grid gave consistent constraints on α, we repeated modeling with a narrowed starting range of α. We show the resulting χ² map in Fig. 3.

We identified three possible topologies, depending on the source trajectory with respect to the position of the caustic (Gaudi & Gould 1997). The inner-large q topology is in the case in which the trajectory crosses between the central and the planetary caustic; outer-large q, in which the trajectory passes behind the planetary caustic with respect to the central caustic; and an intermediate case, small q, in which the trajectory passes through the planetary caustic, possible only in the case of a small mass ratio between the lenses (Fig. 4). The small q topology should suffer the same “inner-outer” degeneracy as in the case of large q topologies. However, the absence of a sudden and strong brightening typical of caustic-crossing in the observed light curve implies that if the trajectory of the source indeed crosses the planetary caustic, the projected size of the source must be comparable to the size of the caustic. Thus, the inner-outer degeneracy in the small q topology does not result in two separated groups of solutions due to the finite-source effect. Finally, we fitted each of the three topologies separately, with all parameters freed. For s and q the starting values were taken from the lowest-χ² models and these parameters were constrained to the ranges marked in Fig. 3.

Fig. 3 Δχ2 map in (log s, log q) parameter space for OGLE-2016-BLG-1609. We named the three visible topologies according to the projected separation: inner-large q (labeled A and delineated by the dashed blue line), small q (B, dot-dashed orange line), and outer-large q (C, dotted green line). Crosses mark the lowest-χ2 models for each topology, and boxes represent the s and q limits used in modeling.

Fig. 4 Source trajectory and caustics of OGLE-2015-BLG-1609 for different binary lens topologies: inner-large q (blue line), small q (orange line), and outer-large q (green line). The central caustics of the small q and inner-large q models at the point (θx/θE, θy/θE) = (0, 0) are approximately point-like. For clarity, we plot only the small q topology trajectory, as the other trajectories differed only slightly. The orange circle represents the size of the source in the case of the small q topology.

3.2 2L1S model with the parallax effect

To obtain an estimate of the absolute values of the lens masses, additional effects in the light curve have to be considered. Typically, the most prominent effect is the microlensing parallax. It is a deviation from the straight trajectory of the source in the lens plane caused by the Earth’s orbital motion. This effect can be used as a source of information on the relative movement of the source and the lens. We introduced this effect into our models by adding north and east components of the microlensing parallax vector, π_E,N and π_E,E, as parameters.

Because the parallax effect breaks the symmetry between positive u₀ and negative u₀ topologies, we additionally split our models between positive and negative u₀ values. The introduction of the parallax effect lowered the χ² statistic of our models. However, for all considered topologies, the obtained absolute values of the north component of the parallax vector were exceptionally large (π_E,N = $-{1.85}_{-0.45}^{+0.61}$, for the inner-large q topology with positive u₀), as shown in Fig. 5. These results contrast with the predictions from the Galactic models, which suggest |π_E,N| ⪅ 0.5 (see Fig. 6). Additionally, fitting the models separately to OGLE and MOA datasets resulted in distinctly different π_E posteriors (Fig. 5).

3.3 Galactic prior

For the new OGLE reductions, we repeated additional analyses conducted for the initial reductions (including the addition of the xallarap effect, the 2L2S models, and the Gaussian-process-correlated noise modeling). However, none of these approaches resulted in reliable improvements in χ² or in the values of π_E For these reasons, we incorporated priors on π_E,N and π_E,E in our modeling, based on the Galactic model from Koshimoto et al. (2021a) and Koshimoto & Ranc (2022). We estimated the input event parameters for the Galactic model by summing the distributions of all three 2L1S u₀+ topologies (inner-large q, small q, outer-large q), and considered the 3σ posterior scatter of the selected parameters, as summarized in Table 3. With those settings, we simulated a sample of 5 × 10⁷ microlensing events with the genulens code (Koshimoto & Ranc 2022). As priors for the Bayesian inference we used the resulting samples of π_E,N and π_E,E, which were broadened in order to obtain proper sampling of wings of the posterior distribution. We achieved that by averaging simulated samples with normal distributions with μ = 0.05, σ = 0.30 and μ = 0.03, σ = 0.28, respectively. These values of μ and σ were determined by fitting Gaussian functions to the simulated π_E,N and π_E,E distributions, with the resulting σ values scaled by a factor of 4. Figure 6 presents the priors alongside the original probability distribution functions of the simulated events. The parameters resulting from the modeling are summarized in Table 4, and the obtained distributions of π_E values are included in Fig. 5.

3.4 Binary-source model

The anomaly in the OGLE-2015-BLG-1609 light curve lacks characteristic features associated with caustic crossing that could unambiguously determine its planetary nature. Therefore, this brightening can alternatively be explained by the presence of a secondary star in the source system. In the single-lens binary-source (1L2S) model, effective magnification is given by the normalized sum of the magnifications of the two sources: $$A=\frac{A_1{f}_{\mathrm{S},1}+A_2{f}_{\mathrm{S},2}}{f_{\mathrm{S},1}+f_{\mathrm{S},2}},$$(6)

where A_i denotes the lensing magnification of each source star with a flux F_i. A_i is given analogously to Eqs. (3) and (4), under the assumption that the transverse speed of the two sources with respect to the lens is the same. Thus, the inclusion of a secondary source doubles the parameters associated with the source and its trajectory (i.e., u_0,1, u_0,2, t_0,1, t_0,2, ρ₁, ρ₂, f_S,1, and f_S,2). The 1L2S model suffers from an analogous degeneracy for positive and negative values of the impact parameter, similar to the single-source interpretation, resulting in two degenerate solutions: (u_0,1+, u_0,2+) and (u_0,1+, u_0,2−), where the signs indicate positive and negative values. The inclusion of the parallax effect additionally breaks the symmetry, leading to four possible solutions: (u_0,1+, u_0,2+),(u_0,1+, u_0,2−), (u^0,1+, u_0,2+), and (u_0,1+, u_0,2−).

Initially, as for binary-lens model, we fitted the binary-source model without accounting of parallax effect separately to the MOA and OGLE data. This model fitted both datasets fairly well, having a significantly lower-χ² statistic than the 1L1S model but still slightly higher than the binary-lens interpretation, in the case of the OGLE data, as presented in Table 2. The MOA dataset is missing observations from two nights at the beginning of the anomaly, leading to weaker constraints on the 1L2S model and a slightly lower-χ² statistic compared to the 2L1S model. Nevertheless, in all cases, the recovered parameters of the models were difficult to justify astrophysically. The timescales of magnification for the two sources are drastically different. For the primary source, it corresponds to the timescale of the entire event, while for the secondary source, it is the timescale of the anomaly. Since t_E is common for both sources, reproducing this difference requires the secondary source to pass significantly closer to the lens than the primary source (|u_0,2| ≈ 0.007) and to be substantially dimmer than the primary source. The absence of a sharp brightening during the anomaly in the light curve further suggests a strong influence of finite-source effects, placing tight constraints on the diameter of the secondary source. In conclusion, in the 1L2S model, the secondary source must be dimmer but not necessarily smaller than the primary source. For the (u_0,1+, u_0,2−) model fitted to the OGLE dataset, ratio of the angular diameters of the sources ρ₁/ρ₂ = ${1.14}_{-0.79}^{+0.85}$ and their brightens in the I-band f_S,1,I/ f_S,2,I = ${195}_{-19}^{+17}$, for MOA dataset these values are ρ₁/ρ₂ = ${2.4}_{-1.5}^{+2.7}$ and f_S,1,R/ f_S,2,R = ${288}_{-51}^{+53}$. Taking into account that the photometric measurements of the source reported in the internal OGLE maps are I = 16.558 ± 0.012 mag and (V − I) = 2.17 ± 0.25 mag, one can conclude that the primary, brighter source is a Bulge giant star. Therefore, the secondary source cannot be simultaneously approximately 200 times dimmer and of comparable size that the primary source. To address this astrophysical inconsistency in the 1L2S models, we assumed that the two sources have similar effective temperatures and included an additional Gaussian prior in our sampling, constraining the ratio of sources flux to the ratio of their diameters f_S,1/ f_S,2= 𝒩((ρ₁/ρ₂)², 10). This resulted in models that are astrophysically justifiable but significantly worsened the χ² statistic of the fits (Table 2). Finally, we attempted to fit the 1L2S model with the parallax effect simultaneously to all available datasets (collected by the surveys and follow-up projects). However, this neither resolved the astrophysical inconsistency in the secondary source properties nor improved the χ² statistic in the sampling that included the mentioned size-diameter prior.

For all the reasons described above, we decided not to consider the binary-source explanation in further analysis.

Fig. 5 Comparison of the posterior distribution of the microlensing parallax vector, πE, for an inner-large q, u0+ topology fitted to different combinations of datasets (from the left): only OGLE data, only MOA data, combined OGLE and MOA datasets, and OGLE data alone with an incorporated prior on πE based on a Galactic model. The orange region marks the 1σ contour of the two-dimensional distribution, and the turquoise region marks the 2σ contour.

Fig. 6 Probability density functions of the two components of the microlensing parallax vector (πE,N in the left panel and πE,E in the right panel), calculated from the Galactic model (Koshimoto et al. 2021a). The dot-dashed turquoise line represents the original function from the simulated events, while the dot-dashed orange function is broadened. The broadened functions are used as priors in the Bayesian inference.

4 Final multi-instrument data modeling: OGLE, MOA, RoboNet, and MiNDSTEp

4.1 Microlensing parameters

To conclude the analysis of OGLE-2015-BLG-1609 and obtain trustworthy parameters of the system, we fitted 2L1S models with the parallax effect simultaneously to all available data from the four mentioned projects: OGLE, MOA, RoboNet, and MiNDSTEp. In all three topologies (inner-large q, small q, and outer-large q), the models with positive values of u₀ had lower values of χ² statistics. These three models are plotted in Fig. 1, and Fig. 2 shows the same models zoomed in on the planetary anomaly. The corresponding parameters of these models, along with their χ² statistics, are listed in Table 5. Error bars scaling parameters for corresponding topologies are summarized in Table 6.

4.2 Physical parameters

4.2.1 Source star

Based on the results from the microlensing modeling, the source star is most likely a typical Galactic bulge red giant. In contrast, the blend is a bluer main-sequence star (see Fig. 7). To acquire extinction parameters, we compared the position of the red clump (RC) giants centroid on the CMD to its intrinsic color and brightness. Based on the OGLE-IV photometric maps for the Galactic bulge, calibrated to the standard systems, we selected stars within 2′ of the location of the event and lying on the red giant branch (i.e., with (V − I) > 2.4 mag and 16.1 < I < 16.7 mag; see Fig. 7). This criterion provided a sample of N_stars = 543. Using these stars, we fitted the luminosity function parameterized as in Nataf et al. (2013): $$\begin{array}{rl}N(I)dI=& A\exp \left[B(I-I_{\mathrm{R}\mathrm{C}})\right]\\ & +\frac{N_{\mathrm{R}\mathrm{C}}}{\sqrt{2\pi }{\sigma }_{\mathrm{R}\mathrm{C}}}\exp [-\frac{{(I-I_{\mathrm{R}\mathrm{C}})}^2}{2{\sigma }_{\mathrm{R}\mathrm{C}}^2}]\\ & +\frac{N_{\mathrm{R}\mathrm{G}\mathrm{B}\mathrm{B}}}{\sqrt{2\pi }{\sigma }_{\mathrm{R}\mathrm{G}\mathrm{B}\mathrm{B}}}\exp [-\frac{{(I-I_{\mathrm{R}\mathrm{G}\mathrm{B}\mathrm{B}})}^2}{2{\sigma }_{\mathrm{R}\mathrm{G}\mathrm{B}\mathrm{B}}^2}]\\ & +\frac{N_{\mathrm{A}\mathrm{G}\mathrm{B}\mathrm{B}}}{\sqrt{2\pi }{\sigma }_{\mathrm{A}\mathrm{G}\mathrm{B}\mathrm{B}}}\exp [-\frac{{(I-I_{\mathrm{A}\mathrm{G}\mathrm{B}\mathrm{B}})}^2}{2{\sigma }_{\mathrm{A}\mathrm{G}\mathrm{B}\mathrm{B}}^2}],\end{array}$$(7)

$$N_{\mathrm{R}\mathrm{G}\mathrm{B}\mathrm{B}}=0.201\times N_{\mathrm{R}\mathrm{C}},$$(8) $$N_{\mathrm{A}\mathrm{G}\mathrm{B}\mathrm{B}}=0.028\times N_{\mathrm{R}\mathrm{C}},$$(9) $$I_{\mathrm{R}\mathrm{G}\mathrm{B}\mathrm{B}}=I_{\mathrm{R}\mathrm{C}}+0.737,$$(10) $$I_{\mathrm{A}\mathrm{G}\mathrm{B}\mathrm{B}}=I_{\mathrm{R}\mathrm{C}}-1.07,$$(11) $${\sigma }_{\mathrm{R}\mathrm{G}\mathrm{B}\mathrm{B}}={\sigma }_{\mathrm{A}\mathrm{G}\mathrm{B}\mathrm{B}}={\sigma }_{\mathrm{R}\mathrm{C}},$$(12) where I, σ, and N are the mean magnitude, magnitude dispersion, and number of stars, respectively. The subscripts AGBB denote the asymptotic giant branch bump and RGBB the red giant branch bump. Following Nataf et al. (2013, 2016), we imposed priors summarized in Table 7.

To sample the posterior distribution, we used the emcee code. Assuming a Poisson distribution of observed stars, we employed the likelihood function in the form $$\ln \mathcal{L}\propto -N_{\mathrm{e}\mathrm{x}\mathrm{p}}+\sum _i\ln N\left(I_i\right)dI.$$(13)

The fitting resulted in I_RC = 16.405 ± 0.47 mag and (I − V)_RC = 2.93 ± 0.25 mag. As seen in the color-magnitude diagram (Fig. 7), interstellar extinction within the selected 2′ radius is highly differential, which is typical for the Galactic bulge region (Nataf et al. 2013). We estimated the reddening parameters using the values of intrinsic apparent magnitude of the RC calculated by Nataf et al. (2016), $$I_{\mathrm{R}\mathrm{C},0}=14.3955-0.0239l+0.0122|b|=14.3350\mathrm{m}\mathrm{a}\mathrm{g}\text{\hspace{0.17em}and\hspace{0.17em}}$$(14) $$(V-I{)}_{\mathrm{R}\mathrm{C},0}=1.06\mathrm{m}\mathrm{a}\mathrm{g},$$(15)

and found $$A_{\mathrm{I}}=2.07\pm 0.47\mathrm{m}\mathrm{a}\mathrm{g}\text{\hspace{0.17em}and\hspace{0.17em}}E(V-I)=1.87\pm 0.25\mathrm{m}\mathrm{a}\mathrm{g}.$$(16)

We combined the reddening with the posterior distribution of the source flux from the microlensing modeling. From this, we obtained the de-reddened color and brightness of the source star. For the inner-large q, u₀+ topology we find $$I_0=15.211\pm 0.058\mathrm{m}\mathrm{a}\mathrm{g}\text{\hspace{0.17em}and\hspace{0.17em}}(V-I{)}_0=0.96\pm 0.25\mathrm{m}\mathrm{a}\mathrm{g}.$$(17)

Fig. 7 Color–magnitude diagram for stars in the OGLE-IV data within 2′ of the microlensed star in the OGLE-2015-BLG-1609 event. The red circle marks the position of the RC centroid, while the turquoise and purple circles mark the positions of the source and blend, respectively.

4.2.2 Source radius

As the color–angular size relations are typically better constrained for the K-band filter, we transformed (V − I)₀ of the source to (V − K)₀ using color-color relations from Bessell & Brett (1988). To calculate the source angular radius corrected for the limb-darkening, we used relations from Adams et al. (2018): $$\log {\theta }_{\ast ,\mathrm{L}\mathrm{D}}=0.562\pm 0.009+0.051\pm 0.003(V-K{)}_0-0.2K_0-\log 2.$$(18)

We assumed that in this relation, as well as in the following ones Eqs. ((20), (21)), the posterior distribution of a given factor is approximated as a normal distribution. Taking the calculated distribution of I₀ and (V − I)₀ for the inner-large q, u₀+ model provided $${\theta }_{\ast ,\mathrm{L}\mathrm{D}}={3.84}_{-0.56}^{+0.85}\mu \mathrm{a}\mathrm{s}.$$(19)

4.2.3 Limb-darkening coefficients

To find the effective temperature of the source star, we used the empirical color relations from Houdashelt et al. (2000) for cold giant stars: $$\begin{array}{rl}{T}_{\text{eff\hspace{0.17em}}}& =8556.22\pm 204.27-5235.57\pm 352.83(V-I{)}_0\\ & +1471.09\pm 148.20(V-I{)}_0^2={4957}_{-720}^{+744}\mathrm{K}.\end{array}$$(20)

Estimation of surface gravity was obtained from the Berdyugina & Savanov (1994) relation for G-K giants and subgiants, assuming a solar metallicity for the source: $$\begin{array}{rl}\log g=& 8.0\pm 0.04T_{\text{eff\hspace{0.17em}}}+0.31\pm 0.04M_{0,\mathrm{V}}\\ & +0.27\pm 0.11[\mathrm{F}\mathrm{e}/\mathrm{H}]-27.15\pm 2.1=3.1\pm 2.2,\end{array}$$(21)

where we used M_0,V = 1.57 ± 0.27 mag. To derive the limb-darkening coefficients, we used tables from Claret & Bloemen (2011) and assumed solar metallicity and turbulent velocity of 2 km s⁻¹. Those parameters have minor impact on the values of coefficients. For bands in which observational data were collected, we obtained $$u_V=0.9198,u_R=0.6535,u_I=0.5597,u_{z^{\mathrm{\prime }}}=0.5131.$$(22)

4.2.4 Lens masses and distance

We obtained the angular Einstein radius based on the scaled source radius from finite-source effect modeling and calculated the observed source radius: $${\theta }_{\mathrm{E}}=\frac{{\theta }_{\ast ,\mathrm{L}\mathrm{D}}}{\rho }={0.46}_{-0.24}^{+1.06}\mathrm{m}\mathrm{a}\mathrm{s}.$$(23)

Combining this with constraints on the microlensing parallax led to the mass measurement of the host and the companion (Gould 2000): $$M_{\mathrm{h}}=\frac{{\theta }_{\mathrm{E}}}{\kappa {\pi }_{\mathrm{E}}}={0.17}_{-0.12}^{+0.63}{\mathrm{M}}_{\odot }\text{\hspace{0.17em}and\hspace{0.17em}}{M}_{\mathrm{c}}=q{M}_{\mathrm{h}}={0.24}_{-0.17}^{+0.90}{\mathrm{M}}_{\mathrm{J}},$$(24)

where κ ≡ 4G/(c²au) ≃ 8.144 mas ${\mathrm{M}}_{\odot }^{-1}$. Figure 8 presents the full distributions of ρ from the MCMC samples and the corresponding distributions of M_h. Additionally, the source-lens relative parallax is given by $${\pi }_{\mathrm{r}\mathrm{e}\mathrm{l}}={\pi }_{\mathrm{E}}{\theta }_{\mathrm{E}}={0.18}_{-0.13}^{+0.51}\mathrm{m}\mathrm{a}\mathrm{s},$$(25)

and the relative proper motion by $${\mu }_{\mathrm{r}\mathrm{e}\mathrm{l}}=\frac{{\theta }_{\mathrm{E}}}{t_{\mathrm{E}}}={5.7}_{-3.0}^{+13.1}\mathrm{m}\mathrm{a}\mathrm{s}{\mathrm{y}\mathrm{r}}^{-1}.$$(26)

Assuming that the source star is at the typical distance for a Galactic bulge star D_S = 8.54 kpc, we estimated the lens distance as $$D_{\mathrm{L}}=\frac{\mathrm{a}\mathrm{u}}{{\pi }_{\mathrm{r}\mathrm{e}\mathrm{l}}+{\pi }_{\mathrm{S}}}={3.3}_{-2.1}^{+2.6}\mathrm{k}\mathrm{p}\mathrm{c},$$(27)

where π_S = 1 / D_S. All the values in this section were calculated for the inner-large q, u₀+ topology. Results for all the considered models are summarized in Table 8.

Fig. 8 Distributions of the radius of the source (ρ) from MCMC samples and calculated distributions of the mass of the host (Mh) of OGLE-2015-BLG-1609 for different binary lens topologies: inner-large q (blue), small q (orange), and outer-large q (green).

4.3 Comparison of solutions

Table 5 includes the χ² values for all three topologies. There is only a minor difference in χ² between the inner-large q and outer-large q topologies. The small q topology has a slightly lower χ² (Δχ² ≈ 10). This difference in χ² is traced to observations from a single night (HJD = 2457267; Fig. 2). On that night, the rising slope of the planetary anomaly in the small q model fits all the OGLE data points, unlike the other two topologies and unlike other nights, where similar variations occurred between measurements. Specifically, in the case of the small q topology, the four OGLE observations taken that night contribute χ²/d.o.f. = 0.9/4 to the overall χ²/d.o.f. = 3239.4/3405. The OGLE observing logs indicate no anomalous conditions during that night, although such a perfect model fit to the data points seems unlikely and indicates overfitting of the small q topology. Unfortunately, no observations from other projects were collected that night to verify the shape of the light curve seen in the OGLE data.

To quantitatively verify the model topologies, we estimated the Bayesian evidence for each of them. This was done using previously simulated microlensing events from the Galactic model of Koshimoto et al. (2021a) with event input properties given in Table 3. Additionally, from the posterior distribution of θ_E for each topologies, we derived the corresponding probability density functions. The evidence was then calculated as $$\mathcal{Z}=\sum _i{p}_i{w}_i\text{\hspace{0.17em}with\hspace{0.17em}}{p}_i=f\left({\theta }_{\mathrm{E},i}\mid \text{model\hspace{0.17em}}\right),$$(28)

where index i denotes i-th simulated event, w_i is its weight, and f is derived probability density function for θ_E for a given model. Evidence values are listed in Table 8. Based on these estimates, only the small q topologies can be considered less probable, while the differences between the outer-large q and inner-large q topologies remains inconclusive.

Considering that 80 M_J is the mass limit for sustaining hydrogen fusion into helium in a stellar core, all topologies suggest that the lens could be a substellar, brown dwarf object. For the two topologies with higher evidence, outer-large q and inner-large q, this probability is estimated to be 49% and 25%, respectively, as shown in Table 8. Using Bayesian evidence values as proportionality indicators, we estimated that the overall probability of the lens being under the 80 M_J mass limit is 34% in the case of the u₀+ models. The small q topology provides the best fit to the data and imposes the strongest constraints on ρ, which in turn constrains the host mass (Fig. 8), yielding an extraordinary host mass of M_h = ${0.025}_{-0.012}^{+0.050}$ M_⊙.

OGLE-2015-BLG-1609L is the 16th known system found through any of the three major planetary detection methods (transit, radial velocity, and microlensing) whose host is suspected to fall within the brown dwarf mass range (Han et al. 2024). All other systems were also detected through microlensing, which results from the unique capabilities of this method in detecting faint objects. Furthermore, the low mass of the lens suggests that the blended flux cannot be fully explained by the lens’s light contribution.

4.4 Prospects for follow-up observations

To better characterize the lens system and more accurately constrain the lens mass, additional high-angular-resolution observations will be needed. Out of the three solutions, only the small q solution has well constrained ρ and, consequently, μ_rel. However, in this solution the μ_rel is very small and the lens is furthest away and of the lowest mass. The combination of the values of these parameters make the lens undetectable in foreseeable future. On the other hand, both the inner-large q and outer-large q solutions have larger μ_rel and closer and more massive lenses. We predict that in 2027 the lens-source separation will be ${68}_{-36}^{+157}$ mas for the inner-large q solution. This separation makes the lens potentially detectable with dedicated adaptive optics instruments, which have the angular resolution of of 60 mas (Vandorou et al. 2023)

Alternatively, use of dual-field interferometer instrument GRAVITY is not possible because of a lack of a bright star (K ≲ 10 mag) for fringe tracking within 30″ of the event’s sky location (Mróz et al. 2025).

5 Summary

We analyzed the microlensing event OGLE-2015-BLG-1609, which shows a planetary anomaly prominent in the surveys data alone. We used a series of modeling approaches, including imposing Galactic model priors, to mitigate challenges such as the unconstrained parallax in the survey data. After incorporating data from the MOA and follow-up projects and refining the OGLE reductions, we identified three possible model topologies for the microlensing system. The two most likely ones have indistinguishable Bayesian evidence. Ultimately, our analysis points to a 34% probability that the lens system is a planet-hosting brown dwarf. We conclude that with a relative proper motion on the order of μ_rel = 5 mas⁻¹, if the lens’s flux is not resolved in the near future by high-angular-resolution instruments, the probability of the lens being a substellar object will increase.

Acknowledgements

We would like to express our sincere gratitude to the anonymous referee for their careful reading of our manuscript and their insightful comments and suggestions. OGLE Team acknowledges former members of the team, for their contribution to the collection of the OGLE photometric data over the past years. The work by Radosław Poleski was partly supported by the Polish National Agency for Academic Exchange grant “Polish Returns 2019”. The MOA project is supported by JSPS KAKENHI Grant Number JP24253004, JP26247023,JP16H06287 and JP22H00153. YT acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (TS 356/3-1). R.F.J. acknowledges support for this project provided by ANID’s Millennium Science Initiative through grant ICN12_009, awarded to the Millennium Institute of Astrophysics (MAS), and by ANID’s Basal project FB210003. L.M. acknowledges financial contribution from PRIN MUR 2022 project 2022J4H55R.

References

All Tables

All Figures

	Fig. 1 Light curves of the microlensing event OGLE-2015-BLG-1609 with three microlensing model topologies for positive values of the impact parameter, u0. The difference between the topologies is visible only during the planetary anomaly (HJD ~ 2457270; see Fig. 2). The photometry comes from the OGLE, MOA, RoboNet, and MiNDSTEp projects.
In the text

	Fig. 2 Same as Fig. 1 but for the light curve of the planetary anomaly in OGLE-2015-BLG-1609.
In the text

Fig. 3 Δχ2 map in (log s, log q) parameter space for OGLE-2016-BLG-1609. We named the three visible topologies according to the projected separation: inner-large q (labeled A and delineated by the dashed blue line), small q (B, dot-dashed orange line), and outer-large q (C, dotted green line). Crosses mark the lowest-χ2 models for each topology, and boxes represent the s and q limits used in modeling.

In the text

Fig. 4 Source trajectory and caustics of OGLE-2015-BLG-1609 for different binary lens topologies: inner-large q (blue line), small q (orange line), and outer-large q (green line). The central caustics of the small q and inner-large q models at the point (θx/θE, θy/θE) = (0, 0) are approximately point-like. For clarity, we plot only the small q topology trajectory, as the other trajectories differed only slightly. The orange circle represents the size of the source in the case of the small q topology.

In the text

Fig. 5 Comparison of the posterior distribution of the microlensing parallax vector, πE, for an inner-large q, u0+ topology fitted to different combinations of datasets (from the left): only OGLE data, only MOA data, combined OGLE and MOA datasets, and OGLE data alone with an incorporated prior on πE based on a Galactic model. The orange region marks the 1σ contour of the two-dimensional distribution, and the turquoise region marks the 2σ contour.

In the text

Fig. 6 Probability density functions of the two components of the microlensing parallax vector (πE,N in the left panel and πE,E in the right panel), calculated from the Galactic model (Koshimoto et al. 2021a). The dot-dashed turquoise line represents the original function from the simulated events, while the dot-dashed orange function is broadened. The broadened functions are used as priors in the Bayesian inference.

In the text

	Fig. 7 Color–magnitude diagram for stars in the OGLE-IV data within 2′ of the microlensed star in the OGLE-2015-BLG-1609 event. The red circle marks the position of the RC centroid, while the turquoise and purple circles mark the positions of the source and blend, respectively.
In the text

	Fig. 8 Distributions of the radius of the source (ρ) from MCMC samples and calculated distributions of the mass of the host (Mh) of OGLE-2015-BLG-1609 for different binary lens topologies: inner-large q (blue), small q (orange), and outer-large q (green).
In the text