MOA-2020-BLG-135Lb: A New Neptune-class Planet for the Extended MOA-II Exoplanet Microlens Statistical Analysis

Gravitational microlensing (Mao & Paczynski 1991) has been solidified as one of the main techniques for detecting planets, being most sensitive to low-mass planets (Bennett & Rhie 1996) that orbit at moderate to large distances from their host star (Gould & Loeb 1992), typically from 0.5 to 10 au, complementing other exoplanet detection methods (Bennett 2008; Gaudi 2012; Batista 2018; Guerrero et al. 2021). The first planetary microlensing event was discovered in 2004 (Bond et al. 2004), and since then, more than 120 exoplanets have been discovered by the method of gravitational microlensing.

The state-of-the-art statistical analysis of planetary signals discovered using gravitational microlensing, Suzuki et al. (2016), implied that cold Neptunes were likely to be the most common type of planets beyond the snow line. This inference was done by discovering a break and likely peak in the planet-to-host-star mass-ratio function for a mass ratio q ∼ 10⁻⁴ when studying the MOA-II microlensing events from 2007 to 2012. At the time of that statistical analysis, it was possible to conclude that while the Suzuki et al. (2016) sample generally supported the predictions for the planet distribution from core accretion theory population synthesis models for planets beyond the snow line (Ida & Lin 2004; Mordasini et al. 2009), the existence of this Neptune peak in the sample distribution actually added contradictions. These previous models for the planet distribution predicted the existence of a sub-Saturn mass planet desert, which conflicted with the microlensing observations (Suzuki et al. 2016, 2018). Only this year Ali-Dib et al. (2022) published their investigation into the origins of cold sub-Saturns, which concluded that these exoplanets may be more common than what was previously predicted. The recent study by Zang et al. (2022) obtained a uniform distribution in $\mathrm{log}q$, which also indicates the absence of a planetary desert of sub-Saturn mass. It is important to notice that the Suzuki et al. (2016) sample had 30 exoplanets, while Zang et al. (2022) had only 13 planets. The size of the samples makes their expansion crucial for future investigations and, consequently, for a better understanding of the distribution of exoplanets.

In this paper, we present the analysis of the gravitational microlensing event MOA-2020-BLG-135 with a short-term planetary lensing signal in its light curve. This analysis leads us to the discovery of MOA-2020-BLG-135Lb, a new planet detected by MOA-II. This new exoplanet qualifies to be included in the upcoming statistical analysis of cold exoplanets detected by the MOA-II survey, which is the expansion of the Suzuki et al. (2016) sample analysis. This paper, presenting a complete study of this event, is organized as follows. First, we describe the observation of the event in Section 2, and how we obtain our best-fit models and explore the full parameter space in Section 3. Then, we explain the photometric calibration and how we retrieve the source size in Section 4, and present our methodology to estimate the lens's physical properties in Section 5, and discuss them comparing with previous literature in Section 6. Finally, we conclude summarizing our results in Section 7.

The microlensing event MOA-2020-BLG-135 was discovered by the Microlensing Observations in Astrophysics (MOA) Collaboration and first alerted on 2020 July 7. This event was located at the J2000 equatorial coordinates $\left({\rm{R}}.{\rm{A}}.,\ \mathrm{decl}.\right)\,=({17}^{{\rm{h}}}\,{53}^{{\rm{m}}}\,41\mathop{.}\limits^{{\rm{s}}}64,-29^\circ 48^{\prime} \,27\mathop{.}\limits^{^{\prime\prime} }24)$, and Galactic coordinates (l, b) = (0.15598°, -1.95678°) in the MOA-II field "gb5" (Sumi et al. 2013). The MOA observations were performed using the purpose-built 1.8 m wide-field MOA telescope located at Mount John Observatory, New Zealand, and the observations of the field "gb5" were taken with a 15 minute cadence using the MOA-Red filter. The MOA-Red filter corresponds to a customized wide-band similar to a sum of the Kron-Cousins R and I bands, from 600 to 900 nm. Occasional observations from the MOA group were made in the visual band using the MOA-V filter. The photometry in these filters was initially performed in real time by the MOA pipeline (Bond et al. 2001), based on the difference imaging method (Tomaney & Crotts 1996). The data used in this paper are from a rereduction done using the Bond et al. (2017) method, which performs a detrending process to correct for systematic errors and removes correlations in the data that may be present due to variations in the seeing and effects of differential refraction (Bennett et al. 2012; Bond et al. 2017). The Bond et al. (2017) method also provides photometry calibrated to the Optical Gravitational Lensing Experiment phase III project (OGLE-III; Szymański et al. 2011).

The Canada–France–Hawaii Telescope (CFHT), located near the summit of Maunakea in Hawaii, United States, also observed the event in the Sloan Digital Sky Survey (Fukugita et al. 1996) i-band filter. From 2020 March to July, CFHT provided one to two supplementary observations per night toward the Korea Microlensing Telescope Network (KMTNet) high-cadence fields and follow-up observations for high-magnification events (Zang et al. 2021). The CFHT data contributes to the establishment of the baseline brightness of the source after the anomaly, and covers a 2 day gap with a data point at HJD = 2, 459, 041.9. The CFHT data were reduced by a custom difference imaging analysis pipeline (Zang et al. 2018) based on the ISIS package (Alard & Lupton 1998; Alard 2000).

The MOA-2020-BLG-135 event was also alerted one day later by the KMTNet group as KMT-2020-BLG-0579. The KMTNet Collaboration monitored this event with the 1.6 m telescope located at the Siding Spring Observatory, Australia (Kim et al. 2016). Unfortunately, in addition to the KMTNet data not covering the anomaly, evidence of systematics was found in their data. Since the data would not improve the characterization of the main event and the planet, the KMTNet team suggested removing their data from the paper.

As a result of observatories' shutdowns due to the COVID-19 pandemic, data from KMT Cerro Tololo Inter-American Observatory, in Chile, KMT South African Astronomical Observatory, in South Africa, and OGLE, in Chile, could not be taken. These observatories were closed when the microlensing event happened.

Figure 1 shows the three data sets used for the analysis of this event. MOA-Red data and the MOA-V data are displayed respectively in brown and violet colors, and the CFHT-i data are in blue.

Zoom In Zoom Out Reset image size

The light curve for the MOA-2020-BLG-135 event (see Figure 1) looks similar to a Paczyński curve (Paczynski 1986), except for the anomaly in the interval $\mathrm{HJD}^{\prime} \,=[9040.7,9041.5]$ ²¹ observed by both MOA-Red and MOA-V (zoomed in Figure 1). The Paczyński curve assumes a model in which the lens consists of a single star and the radiant flux comes from a single source. We display this curve as a point-source point-lens (PSPL) model in Figure 1. This deviation indicates that the lens may be composed of two masses, in which the less massive lens component can be a planet-mass object. We call this a binary or planetary lens system, depending on the mass ratio, as there are two objects contributing gravitationally as lenses. In Section 3.1, we search for a lensing model explaining the three data sets presented in Section 2 (MOA-Red, MOA-V, and CFHT-i) by exploring the parameter space of possible binary lenses with a single source (2L1S) using the method described in Bennett (2010). In Section 3.2, we discuss the lack of evidence for a binary-source model—with a single lens (1L2S) or a binary lens (2L2S)—after investigating the light curve using the method described in Bennett et al. (2018).

3.1. Single-source Scenario

The Bennett (2010) process uses the image-centered, ray-shooting method (Bennett & Rhie 1996) combined with a custom version of the Metropolis algorithm (Metropolis et al. 1953), which yields a rapid convergence to a χ² minima.

3.1.1. Fit Parameters

The parameters of our model are the Einstein crossing time (t_E); the time at which the separation of the lens and source reaches the minimum (t₀); the minimum angular separation between the source and lens as seen by the observer (u₀); the separation of the two masses of the binary-lens system during the event (s); the counterclockwise angle between the lens-source relative motion projected onto the sky plane and the binary-lens axis (α); the mass ratio between the secondary lens and the primary lens (q); the source radius crossing time (t_*); the source flux for each instrument i (f_s,i ); and the blend flux per instrument i (f_b,i ).

The parameters t_E, t₀, and u₀ are the common parameters for the single-lens model, while s, α, and q are the additional parameters for a binary-lens system model. Both length parameters, u₀ and s, are normalized by the angular Einstein radius θ_E, defined by

$$\begin{eqnarray}&&{\theta }_{{\rm{E}}}=\sqrt{\displaystyle \frac{4{{GM}}_{L}}{{c}^{2}{D}_{S}}\left(\displaystyle \frac{{D}_{S}}{{D}_{L}}-1\right)},\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, M_L is the mass of the lens system, c is the speed of light, D_S is the observer-source distance, and D_L is the observer-lens distance. The source radius crossing time, t_*, is included in the lensing model to take account of finite source effects,

$$\begin{eqnarray}&&{t}_{* }=\rho \ {t}_{{\rm{E}}}=\displaystyle \frac{{\theta }_{* }}{{\theta }_{{\rm{E}}}}{t}_{{\rm{E}}},\end{eqnarray} \tag{ 2 }$$

where ρ is the source angular radius in Einstein units, and θ_* is the source angular radius.

The other two parameters taken into account are the blend flux f_b,i and the source flux f_s,i . As microlensing events are observed in crowded stellar fields, the source is usually blended with other unlensed stars. For this reason, we consider the blend flux. Since the observed brightness has a linear dependence on the blend flux and the source flux, they are treated differently from the other nonlinear fit parameters, as it follows. For every instrument and each set of the previously cited fit parameters, we can find a total flux that minimizes the χ². The total flux F_i (t) for time t for instrument i can be written as

$$\begin{eqnarray}&&{F}_{i}(t)=A(t,x)\ {f}_{s,i}+{f}_{b,i},\end{eqnarray} \tag{ 3 }$$

where A(t, x) is the magnification of the event at any given time and for any given set of nonlinear parameters x = (t_E, t₀, u₀, s, α, q, t_*), f_s,i is the unlensed source flux in a specific passband i, and f_b,i is its excess flux. For many light curves reduced with difference imaging, the blend flux has an arbitrary normalization. Yet, the Bond et al. (2017) method normalizes the total flux to match the flux of the nearest star-like object in the reference frame, and it is the one used for the MOA data.

For this modeling, we do not consider parallax effects because this is a short (t_E ≈ 17 days ≪ 1 month) and faint event. Additionally, its peak of magnification was reached on 2020 July 10, when the Earths instantaneous acceleration toward the projected position of the Sun projected into the lens plane was close to its minimum. These three factors make it very unlikely to detect asymmetric features in the light-curve tails due to parallax effect. Therefore, we do not attempt a parallax measurement for this event.

3.1.2. Exploring the Full Parameter Space

We start modeling by systematically exploring the full parameter space with the grid-search approach described in Bennett (2010). For the first step, we do an initial condition grid search where t_E, t₀, u₀, and t_* are fixed, while s, α, and q are scanned. This is done so we can select the initial conditions. From there, we use the custom version of the Metropolis algorithm, as implemented by Bennett (2010), with initial positions coming from the 13 local solutions obtained from these initial scans. To ensure that no good model is missed in our analysis, we make sure that different mass ratios in the q ∈ [10⁻⁵, 10⁻¹] range are included in our set of initial conditions. Our customized Metropolis algorithm provides a full fit, in which all the parameters are free. Then, we select the best-fit model by looking for the lowest χ², which indicates a model with q ∼ 10⁻⁴.

Interpreting a planetary lensing event with a planet signal is often subject to a close-wide degeneracy, in which a solution with a certain separation s results in a similar model to another model presenting a solution with the same parameters' values but with a separation given by 1/s. This degeneracy arises from the symmetries in the lens equation (Griest & Safizadeh 1998; Dominik 1999). As discussed in Yee et al. (2021), the degeneracy appears even for resonant caustics, which are far from the s ≫ 1 limit in which the symmetries were derived. Therefore, we carefully cover both the close (s < 1) and the wide (s > 1) solutions in our analysis.

To guarantee the exploration of the full parameter space with the Monte Carlo method, we use our customized version of the Metropolis algorithm for both the close and wide solutions. When running the Markov Chain Monte Carlo (MCMC) algorithm with the parameters of the best-fit models as initial inputs, we notice that the proposal distribution function we choose allows each chain to jump back and forth between the wide and close solutions. This happens because both values for the separation s are close enough that the χ² barrier between them has a size encompassed by the possible jumps. To ensure the optimization of the posterior sampling, we conduct MCMC runs adjusting the variables for our diagonalized covariance matrix and combine the results of the runs. For our analysis, we combine only the MCMC runs that jump between the wide and close solution, and vice versa, and those have chains with a similar best-fit model.

Our best-fit planetary light-curve model (2L1S s > 1) is shown in Figure 1, and its parameters are given in Table 1, being the solution referred to as simply "the best-fit model" in this paper. The result of our best-fit model for 2L1S s < 1 is also displayed in Table 1. The median of the marginalized posterior distributions with the 1σ confidence interval (i.e., 68.3%) is displayed in the same table, together with the range for distributions within the 2σ interval (i.e., 95.5%). Figure 2 shows the posterior distribution together with the 1σ, 2σ, and 3σ (i.e., 99.7%) confidence intervals.

Zoom In Zoom Out Reset image size

Table 1. Best-fit Model Parameters for MOA-2020-BLG-135, and Corresponding Medians from the Posterior Distribution

Parameters	Units	2L1S s < 1	2L1S s > 1	MCMC Medians	2σ Range
tE	days	16.905	16.791	${16.85}_{-0.27}^{+0.28}$	16.30–17.41
t0	$\mathrm{HJD}^{\prime}$	9040.48797	9040.48843	${9040.4879}_{-0.0019}^{+0.0020}$	9040.4839–9040.4920
u0		0.03260	0.03287	0.03275 ± 0.00072	0.03134–0.03421
s		0.95472	1.04509	${0.997}_{-0.054}^{+0.096}$	0.916–1.126
α	radians	2.38728	2.37731	${2.383}_{-0.013}^{+0.021}$	2.360-2.433
q	10−4	1.55314	1.13963	${1.52}_{-0.31}^{+0.39}$	1.01–2.47
t*	days	0.14661	0.14186	${0.145}_{-0.011}^{+0.018}$	0.125–0.183
IS		19.009	19.000	19.004 ± 0.045	18.914–19.095
VS		21.137	21.127	21.132 ± 0.045	21.042–21.223
fit χ2		14819.9	14819.7	⋯	⋯

Download table as:  ASCIITypeset image

In Figure 2, a butterfly-wing shape is visible in the posterior distribution for s. This shape shows that both the close and wide solutions were fully explored during our MCMC runs. This is an effect of our algorithm jumping back and forth between both wide and close solutions. Even though there seems to exist two regions in each butterfly wing when considering the 1σ interval, the χ² difference is not big enough to create a barrier that could separate them into four independent regions. Yang et al. (2022) discuss the existence of pairs of close/wide solutions, which was named the "central-resonant degeneracy." The mentioned interconnected four local minima in Figure 2 could be interpreted as two of those pairs. Figures 3(a) and (b) illustrate the close-separation (s < 1) and the wide-separation (s > 1) topology, respectively. For the s > 1 solution, a resonant caustic solution appears as the best solution, while for the s < 1 solution, a nonresonant caustic solution is slightly better. It is not surprising that the s > 1 and s < 1 regions include resonant and nonresonant caustic solutions because the caustic topology is related to light-curve features only in special cases. The q/10⁻⁴ versus s panel of Figure 2 shows the dividing lines for these caustic topologies. For our data points, the best-fit model for the light curve with the planetary separation solution s > 1 is practically the same as the one for the solution s < 1. For s < 1, the anomaly due to the presence of the planet appears to start at $\mathrm{HJD}^{\prime} \approx 9040.85$, finishing at $\mathrm{HJD}^{\prime} \approx 9041.30$, while for s > 1, it appears to start at $\mathrm{HJD}^{\prime} \approx 9040.82$, finishing at $\mathrm{HJD}^{\prime} \approx 9041.34$. The explanation of these small differences can be found in Figure 3, which indicates the source starts and finishes crossing the caustic at slightly different times. Moreover, we can check the great similarity for all the fitting parameters, together with both s being close to the inverse of each other, indicating the historically called close-wide degeneracy solutions.

Zoom In Zoom Out Reset image size

Examples of close-wide model degeneracies that do not obey the s ↔ 1/s relationship predicted by Dominik (1999) are relatively common (e.g., Bennett et al. 2014; Koshimoto et al. 2017), and An (2021) pointed out that many such cases can be explained similarities of local caustic regions even in situations in which the overall caustic shape may not be degenerate at all. Adding to the degeneracy discussion, Zhang et al. (2022) have recently demonstrated that the close-wide degeneracy s ↔ 1/s relationship is only strictly followed for the singular case of u₀ = 0. For the general case of u₀ > 0, the authors proposed an alternative theory named "the offset degeneracy" that predicts a deviation from s ↔ 1/s. The formalism is shown to be mathematically exact in certain limits (Zhang & Gaudi 2022), which includes caustic-crossing events and resonant events. In this work, the different degenerate solutions do not indicate distinct conclusions for the posterior distribution of the physical properties. Therefore, the solutions have a similar lens physical interpretation.

One might wonder about the proximity of the source and the upper cusp at $\mathrm{HJD}^{\prime} \approx 9041.71$ in the lower panel of Figure 3. This upper cusp does not seem to create an extra perturbation, and it is located far from the source, more than one source radius away at $\mathrm{HJD}^{\prime} \approx 9041.71$. It results in no apparent anomaly in the best-fit model in Figure 1.

3.2. Binary-source Possibility

The source angular radius θ_* is not explicitly obtained from the light-curve models presented in Section 3, yet it can be empirically derived if we know the dereddened magnitude and color of the source (Van Belle 1999; Yoo et al. 2004; Bennett 2010). Toward this aim, we correct the extinction and reddening for the source star by using the red clump giants as a reference.

In order to obtain the dereddened color and corrected magnitude of our source star, we first calibrate the instrumental MOA-II magnitudes, MOA-Red and MOA-V, with the OGLE-III catalog considering the following equations from the standard calibration procedure described in Bond et al. (2017):

$$\begin{eqnarray}\begin{array}{rcl}{I}_{{\rm{O}}3} & = & (28.0983\pm 0.0014)+{R}_{{MOA}}\\ & & -(0.20844\pm 0.00087)\times ({V}_{\mathrm{MOA}}-{R}_{\mathrm{MOA}})\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{V}_{{\rm{O}}3} & = & (28.5038\pm 0.0014)+{V}_{\mathrm{MOA}}\\ & & -(0.10746\pm 0.00088)\times ({V}_{\mathrm{MOA}}-{V}_{\mathrm{MOA}}),\end{array}\end{eqnarray} \tag{ 5 }$$

where R_MOA is the MOA-Red filter and V_MOA is the MOA-V filter. The Bond et al. (2017) calibration is the standard MOA calibration procedure and is done using crossmatched stars within the 500'' × 500'' cameo image from our DOPHOT (Schechter et al. 1993) catalog with stars from the OGLE-III catalog. In the final fit, the rms scatters, which are dominated by thousands of faint stars, are rms_I ≈ 0.07 and rms_V ≈ 0.08. We then obtain the magnitude I_O3 and V_O3 (Gould et al. 2010) for our source star. These are displayed in Table 1 as I_S and V_S respectively. Figure 4 shows the calibrated color and magnitude for our source star, compared to the stars within a 90'' radius limit from our target.

Zoom In Zoom Out Reset image size

For the second step, we measure the extinction and reddening of the stars within a 90'' radius limit from our source star as follows. First, we calculate what are the apparent magnitude and color of the red clump, obtaining I_RCG =(16.09 ± 0.05) and (V − I)_RCG = (2.40 ± 0.05) (represented as the red dot in Figure 4). The expected dereddened magnitude of the red clump at a Galactic longitude l = 0.15598° is I_RCG,0 = (14.44 ± 0.04) (Nataf et al. 2013), and the expected color is (V − I)_RCG,0 = (1.06 ± 0.06) (Bensby et al. 2011). Therefore, we can use the calculated apparent magnitude and color of the red clump in comparison with the expected values to obtain the extinction and the color excess of this region in the sky. The extinction and the color excess are A_I =(1.65 ± 0.06) and E(V − I) = (1.34 ± 0.08). These values are reasonable, and can be compared to the ones given by the OGLE-III tool for querying interstellar extinction toward the Galactic bulge, ²² which were A_I = 1.60 and E(V − I) = 1.35 when using the natural neighbor interpolation option.

With our calculated extinction and the color excess, we obtain the corrected magnitude I_S,0 = (17.350 ± 0.078) and dereddened color ${(V-I)}_{{\rm{S}},0}={0.788}_{-0.096}^{+0.095}$ (in Table 3) of our source star.

Table 3. Source and Lens-source Properties

Parameters	Units	MCMC Medians
Source magnitude IS,0		${17.350}_{-0.078}^{+0.078}$
Source magnitude KS,0		${16.45}_{-0.25}^{+0.22}$
Source color (V − I)S,0		${0.788}_{-0.096}^{+0.095}$
Source angular radius θ*	μas	${1.15}_{-0.12}^{+0.13}$
Einstein radius θE	mas	${0.133}_{-0.018}^{+0.019}$
Lens-source proper motion μrel,G	mas yr−1	${2.88}_{-0.40}^{+0.42}$

Download table as:  ASCIITypeset image

Finally, we determine the source size:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left[\displaystyle \frac{2{\theta }_{* }}{\mathrm{mas}}\right]=0.501414\ +\ 0.419685{(V-I)}_{{\rm{S}},0}\ -\ 0.2{I}_{{\rm{S}},0},\end{eqnarray} \tag{ 6 }$$

following Boyajian et al. (2014) analysis for stars with 3900 K < T_eff < 7000 K, and appearing in Bennett et al. (2017). The source angular radius is ${\theta }_{* }={1.15}_{-0.12}^{+0.13}$ μas.

Determining the source size is necessary to calculate the angular Einstein radius θ_E when using Equation (2). The importance of θ_E itself comes from the fact that the main physical properties of any microlensing event (the lens mass M_L, the distance to the lens D_L, the distance to the source D_S, and the lens-source relative proper motion μ_rel,G) are directly related to it. The measurement of θ_E provides one mass-distance relationship, defined in Equation (1), which can be rearranged into

$$\begin{eqnarray}&&{M}_{{\rm{L}}}=0.1228{M}_{\odot }{\left(\displaystyle \frac{{\theta }_{{\rm{E}}}}{1\ \mathrm{mas}}\right)}^{2}{\left(\displaystyle \frac{{\pi }_{\mathrm{rel}}}{1\ \mathrm{mas}}\right)}^{-1},\mathrm{and}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{D}_{{\rm{L}}}=1\ \mathrm{kpc}{\left(\displaystyle \frac{{\pi }_{\mathrm{rel}}}{1\ \mathrm{mas}}+{\left(\displaystyle \frac{{D}_{{\rm{S}}}}{1\ \mathrm{kpc}}\right)}^{-1}\right)}^{-1},\end{eqnarray} \tag{ 8 }$$

where the lens-source relative parallax π_rel is defined as

$$\begin{eqnarray}&&{\pi }_{\mathrm{rel}}=\displaystyle \frac{1\ \mathrm{au}}{{D}_{{\rm{L}}}}-\displaystyle \frac{1\ \mathrm{au}}{{D}_{{\rm{S}}}}.\end{eqnarray} \tag{ 9 }$$

We use the Einstein crossing time t_E and the source radius crossing time t_* from our light-curve model, combined with our calculated source angular radius θ_*, to obtain the angular Einstein radius ${\theta }_{{\rm{E}}}={0.133}_{-0.018}^{+0.019}\ \mathrm{mas}$. Even though we can define the lens-source proper motion μ_rel,G as a function of θ_E, we use the following formula for θ_* and t_* that avoids an increased uncertainty due to the blending degeneracy:

$$\begin{eqnarray}&&{\mu }_{\mathrm{rel},{\rm{G}}}=\displaystyle \frac{{\theta }_{{\rm{E}}}}{{t}_{{\rm{E}}}}=\displaystyle \frac{{\theta }_{* }}{{t}_{* }}.\end{eqnarray} \tag{ 10 }$$

The lens-source proper motion is ${\mu }_{\mathrm{rel},{\rm{G}}}={2.88}_{-0.40}^{+0.42}\ \mathrm{mas}\ {\mathrm{yr}}^{-1}$.

Aiming to provide more information about the source for future high-angular-resolution follow-up observations, we also estimate the magnitude of the source in the K band by using the color transformations presented in Kenyon & Hartmann (1995). The source magnitude is ${K}_{{\rm{S}},0}={16.45}_{-0.25}^{+0.22}$ without extinction. By adding the extinction A_K = 0.2195 (Nishiyama et al. 2009; Gonzalez et al. 2012), we obtain ${K}_{{\rm{S}}}={16.67}_{-0.25}^{+0.22}$ as the magnitude of the source in the K band.

The calculated source and lens-source properties described in this section are in Table 3.

As parallax and lens brightness measurements are missing for the MOA-2020-BLG-135 event, we cannot uniquely determine the lens mass and its distance. Therefore, to estimate the lens properties, we use the Bennett et al. (2014) Galactic model. The strength of this model is that it can incorporate a prior for the Bayesian analysis when estimating the posterior probability distribution of the host mass. It allows us to use a mass function under the most conventional assumption that all stars have an equal planet-hosting probability, or under the assumption that planets are more likely to orbit around more-massive stars, by setting a prior proportional to M.

Usually, in microlensing papers, the estimations for the lens system properties assume that all the stars have equal probability of hosting a planet, which implies a mass function prior uniform in M. Statistical results on exoplanet populations were also obtained under the same assumption (Cassan et al. 2012). Yet, for this paper, we also consider a second scenario in which the probability of hosting a planet scales in proportion to the host-star mass, dP ∝ M dM. Johnson et al. (2007, 2010) found that, for their radial velocity sample, the planet occurrence increases with the stellar mass at fixed planet mass. This is compatible with Nielsen et al. (2019) direct imaging sample, which showed a strong correlation between planet occurrence rate and host-star mass. Moreover, Bhattacharya et al. (2021) identified that the traditional assumption, which considers that all the stars have equal probability of hosting a planet, is not consistent with many microlensing events that have been revisited with the help of the Keck adaptive optics and had their lens object identified using high-angular-resolution observations. Bhattacharya et al. (2021) pointed out that five of the six events with direct measurement of the separation between the source and the lens stars have found a host star more massive than the median predicted under the most conventional assumption, which is the one with a prior uniform in M. The authors also indicated that there is no publication bias for that Keck sample. Certainly, a more extensive sample to state this with more confidence is needed and, in fact, NASA Keck Key Strategic Mission Support and Hubble Space Telescope observing programs will be directly measuring the mass of more microlens host stars. Although planets are more likely to orbit around more-massive stars, we still decide to consider both mass priors for our Bayesian assumptions, the conventional prior uniform in the stellar mass, and a prior that scales in proportion to the stellar mass.

Results of the Bayesian analysis can also be affected when the probability of having a planet depends on the stellar location in our Galaxy. However, it has recently been shown by Koshimoto et al. (2021) that the dependence of the planet-hosting probability on the Galactocentric distance is not large; thus we do not consider such dependence in this paper.

The Bayesian analysis for the lens properties is done using as input the collection of the MCMC runs (discussed in Section 3.1.2), with their fit parameters, along with our calculated angular source radius and extinction (see Section 4). To obtain the magnitude of the lens not only in the I band, but also in the K band, which is useful for future high-angular-resolution follow-up observations, we use the extinction A_K = 0.2195 (Gonzalez et al. 2012; Nishiyama et al. 2009). The Bennett et al. (2014) Galactic model is used two times, with the two different priors. First, we consider a power-law stellar mass function under the standard assumption that all stars have an equal planet-hosting probability, a function of the form dP ∝ M^β dM with β = 0. Then we consider the same function but under the assumption that planets are more likely to orbit around more-massive stars, so a function with β = 1.

Under the assumption that the probability of hosting a planet is the same for all stars, β = 0 ⇒ dP ∝ dM, we estimate that the lens physical properties and their 1σ (68.3%) interval of confidence are ${m}_{\mathrm{planet}}={11.3}_{-6.9}^{+19.2}\ {M}_{\oplus }$ for the planet mass, ${M}_{\mathrm{host}}={0.23}_{-0.14}^{+0.39}\ {M}_{\odot }$ for the host mass, ${D}_{{\rm{L}}}={7.9}_{-1.0}^{+1.0}\ \mathrm{kpc}$ for the distance to the lens, ${a}_{\perp }={1.11}_{-0.20}^{+0.23}\ \mathrm{au}$ for the project separation, ${a}_{3{\rm{D}}}={1.35}_{-0.32}^{+0.75}\ \mathrm{au}$ for the deprojected separation, and ${I}_{{\rm{L}}}={26.0}_{-2.9}^{+2.4}$ and ${K}_{{\rm{L}}}={22.2}_{-2.3}^{+1.9}$ for the lens magnitude. When assuming the mass function prior proportional to M, β = 1 ⇒ dP ∝ M dM, the planet mass estimation is ${m}_{\mathrm{planet}}\,={25}_{-15}^{+22}\ {M}_{\oplus }$, the host mass is ${M}_{\mathrm{host}}={0.53}_{-0.32}^{+0.42}\ {M}_{\odot }$, the distance to the lens is ${D}_{{\rm{L}}}={8.3}_{-1.0}^{+0.9}\ \mathrm{kpc}$, the projected separation is ${a}_{\perp }={1.17}_{-0.20}^{+0.23}\ \mathrm{au}$, the deprojected separation is ${a}_{3{\rm{D}}}\,={1.42}_{-0.32}^{+0.78}\ \mathrm{au}$, and lens magnitude is ${I}_{{\rm{L}}}={23.8}_{-3.1}^{+2.5}$ and ${K}_{{\rm{L}}}={20.4}_{-2.3}^{+2.0}$. We report the probability distributions considering both priors in Table 4.

Table 4. Lens Physical Properties Derived from the Bennett et al. (2014) Galactic Model

Parameters	Units	Prior Uniform in M, dP ∝ dM	2σ Range (95.5%)	Prior Proportional to M, dP ∝ M dM	2σ Range (95.5%)
Planet mass mplanet	M⊕	${11.3}_{-6.9}^{+19.2}$	2.1-57.3	${25}_{-15}^{+22}$	4–70
Host mass Mhost	M⊙	${0.23}_{-0.14}^{+0.39}$	0.05-1.07	${0.53}_{-0.32}^{+0.42}$	0.08–1.22
Lens distance DL	kpc	${7.9}_{-1.0}^{+1.0}$	5.9-9.7	${8.3}_{-1.0}^{+0.9}$	6.4–9.9
Projected separation a⊥	au	${1.11}_{-0.20}^{+0.23}$	0.74-1.58	${1.17}_{-0.20}^{+0.23}$	0.80–1.64
Deprojected separation a3D	au	${1.35}_{-0.32}^{+0.75}$	0.82-4.72	${1.42}_{-0.32}^{+0.78}$	0.88–4.97
Lens magnitude IL		${26.0}_{-2.9}^{+2.4}$	20.2-36.7	${23.8}_{-3.1}^{+2.5}$	19.5–28.7
Lens magnitude KL		${22.2}_{-2.3}^{+1.9}$	17.6-32.0	${20.4}_{-2.3}^{+2.0}$	17.1–24.2

Download table as:  ASCIITypeset image

Figure 5 shows the probability distribution of the planet and host masses, the distance to the lens system, their projected separation, and the lens magnitudes in both the I band and K band, under the assumption of equal planet-hosting probability, dP ∝ dM. Figure 6 shows the same results, but for the probability scaling in proportion to the host mass, dP ∝ M dM. It is interesting to note that, even though the host mass and the planetary mass medians seem to be almost the double when comparing the results when the prior is uniform in M to the results when the prior is proportional to M, the range for the masses are not that different when considering the 2σ (i.e., 95.5%) confidence interval.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our light-curve analysis for the event MOA-2020-BLG-135 leads to the discovery of MOA-2020-BLG-135Lb, a new Neptune-class planet. This analysis yields a planet–host mass ratio of $q={1.52}_{-0.31}^{+0.39}\times {10}^{-4}$, and separation s ≈ 1. It is important to mention that with each MCMC chain we were able to sample all the modes of the posterior distribution due to the closeness of our close-separation (s < 1) and wide-separation (s > 1) solutions.

In Sections 4 and 5, we determine the source and the lens magnitude in the K band, aiming to anticipate results for high-angular-resolution follow-up observations. The source magnitude with added extinction was computed to be ${K}_{{\rm{S}}}\,={16.67}_{-0.25}^{+0.22}$. Under the assumption that all stars have equal planet-hosting probability, the lens magnitude is expected to be in the range 19.9−24.1 for the 68.3% (1σ) confidence interval, and 17.6–32.0 for the 95.5% (2σ) interval, with median 22.2. One might wonder whether the lens star looks faint in comparison with the source star when considering only the median and the 1σ interval, which can indicate that its observation might be challenging in near future. However, for this assumption, when considering the 2σ interval of confidence, the result looks more encouraging. The event MOA-2007-BLG-400 had its lens object successfully identified using high-angular-resolution observations, and the magnitude for the source and lens in the Keck K band were, respectively, 16.43 ± 0.04 and 18.93 ± 0.08 (Bhattacharya et al. 2021), which are not that different from our event, when considering the brightest possible lens magnitude. Moreover, on the assumption that the probability of hosting a planet scales in proportion to the stellar mass, the scenario gets even better, the lens magnitude is expected to be in the range 18.1–22.4 for the 68.3% (1σ) confidence interval, and 17.1–24.2 for the 95.5% (2σ) interval, with median 20.4. Knowing that the microlensing events tend to have a host star more massive than their median predictions from only the light-curve analysis, direct detection of the lens star looks promising for our event MOA-2020-BLG-135.

Another notable detail from our results is that a planet with a mass ratio $q={1.52}_{-0.31}^{+0.39}\times {10}^{-4}$ seems to be exactly where previous core accretion theories predicted a Neptune desert (Ida & Lin 2004; Mordasini et al. 2009). It also lies exactly in the peak of the planet-to-star mass-ratio distribution measured by the state-of-the-art statistical analysis of planets detected by gravitational microlensing (Suzuki et al. 2016). Recently, Zang et al. (2022) reported the analysis of a statistical planetary sample with 13 homogeneously selected planets observed in 2019 by the KMTNet Collaboration. The paper suggests that the mass-ratio function may not decrease rapidly below the Suzuki et al. (2016) mass-ratio break (q_br ∼ 10⁻⁴), instead, there may be a uniform distribution in $\mathrm{log}q$. Once more, the Suzuki et al. (2016) contained 30 planets and Zang et al. (2022) contained 13 planets. To better understand this planet-to-star mass-ratio distribution, we need a larger microlens exoplanet sample. The MOA Collaboration has been working to obtain this extended sample, and will have more than 50 new planets, including the planet presented in this paper. Meanwhile, the KMTNet Collaboration is also working on a sample, predicting about 120 planets from 2016 to 2019. It is important to notice that the specifications of the telescopes and their cadence and sensitivity are different; therefore, both analyses are necessary.

In Section 5, we discussed the Bayesian priors used for the analysis of the lens system properties subject to the microlensing light-curve constraints. In almost all previous analyses with similar light-curve constraints, it has been assumed that all host stars have an equal probability to host the planet with the measured mass ratio ($q={1.52}_{-0.31}^{+0.39}\times {10}^{-4}$ for this event). However, Bhattacharya et al. (2021) have shown that higher mass host stars appear to be more likely for the planets that the microlensing method is sensitive to. This is somewhat similar to earlier findings from radial velocity (Johnson et al. 2007, 2010) and direct detection (Nielsen et al. 2019) surveys. However, both of these analyses considered planet-hosting probabilities at a fixed planet mass instead of a fixed mass ratio. Since lower mass ratio planets are more common (Suzuki et al. 2016), for mass ratios ≥10⁻⁴, the demonstration that higher mass hosts are more likely for a fixed mass ratio is stronger than the same claim at a fixed planet mass. If the planet-hosting probability was independent of host mass at a fixed mass ratio, the hosting probability would be higher for higher mass hosts at a fixed planet mass, since this implies a larger mass ratio, q, for the lower mass hosts (as long as ≥10⁻⁴).

We have presented the discovery of MOA-2020-BLG-135Lb, a new Neptune-class planet uncovered by the light-curve analysis for the microlensing event MOA-2020-BLG-135. By applying the Bennett (2010) process, our analysis has revealed a planet–host mass ratio of $q={1.52}_{-0.31}^{+0.39}\times {10}^{-4}$, and separation s ≈ 1.

To estimate the lens system properties for MOA-2020-BLG-135, we have conducted a Bayesian analysis using the Bennett et al. (2014) Galactic model. When considering that all stars have equal probability of hosting a planet, using a mass function prior uniform in M, we have found a planet mass of ${m}_{\mathrm{planet}}={11.3}_{-6.9}^{+19.2}\ {M}_{\oplus }$, a host star with a mass of ${M}_{\mathrm{host}}\,={0.23}_{-0.14}^{+0.39}\ {M}_{\odot }$ and magnitude ${K}_{{\rm{L}}}={22.2}_{-2.3}^{+1.9}$, located at a distance ${D}_{{\rm{L}}}={7.9}_{-1.0}^{+1.0}\ \mathrm{kpc}$. Under the assumption that the hosting probability scales in proportion to the stellar mass, we have estimated ${m}_{\mathrm{planet}}={25}_{-15}^{+22}\ {M}_{\oplus }$, ${M}_{\mathrm{host}}={0.53}_{-0.32}^{+0.42}\ {M}_{\odot }$, ${K}_{{\rm{L}}}={20.4}_{-2.3}^{+2.0}$, and ${D}_{{\rm{L}}}={8.3}_{-1.0}^{+0.9}\ \mathrm{kpc}$.

The previous results from radial velocity and direct imaging (Johnson et al. 2007, 2010; Nielsen et al. 2019) considered planet-hosting probabilities for a fixed planet mass, whereas the MOA-II exoplanet microlens statistical analysis considers the mass ratio. MOA-2020-BLG-135Lb is an important detection for completeness of the extended MOA-II exoplanet microlens statistical sample. Additionally, high-angular-resolution follow-up observations for this event are certainly recommended in the future for restricting the mass values of the host star and its planet.

We thank the government of New Zealand for their strict strategy against COVID-19, which gave MOA minimal down time during the pandemic. The observatory was kept open for most of the observing season, allowing us to observe this planetary event in 2020. We must also thank the KMTNet Collaboration for the valuable paper discussion, especially Jennifer Yee, Cheongho Han, and Andrew Gould. We would also like to thank Keming Zhang for his helpful conversation about his two papers.

Work by S.I.S., D.P.B., A.B., A.V., and G.R. was supported by NASA under award number 80GSFC21M0002. The MOA project is supported by JSPS KAK-ENHI grant No. JSPS24253004, JSPS26247023, JSPS23340064, JSPS15H00781, JP16H06287,17H02871, and 19KK0082. Work by C.R. was supported by the Alexander von Humboldt Foundation, and was partly carried out within the framework of the ANR project COLD-WORLDS supported by the French Agence Nationale de la Recherche with the reference ANR-18-CE31-0002. S.M. acknowledges support from JSPS KAKENHI grant No. 21J11296. This research uses data obtained through the Telescope Access Program (TAP), which has been funded by the TAP member institutes. W.Zang, S.M. and W.Zhu acknowledge support by the National Science Foundation of China (grant No. 12133005). W.Zhu acknowledge the science research grants from the China Manned Space Project with No. CMS-CSST-2021-A11. This work is partly based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada–France–Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the French Centre National de la Recherche Scientifique (CNRS), and the University of Hawaii. This research is supported by Tsinghua University Initiative Scientific Research Program (Program ID 2019Z07L02017).

Software: Matplotlib (Hunter 2007), MOAna (Ranc 2020), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020).