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

We report the light-curve analysis for the event MOA-2020-BLG-135, which leads to the discovery of a new Neptune-class planet, MOA-2020-BLG-135Lb. With a derived mass ratio of q


INTRODUCTION
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-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 most recent statistical analysis of planetary signals discovered using gravitational microlensing implied that cold Neptunes were likely to be the most common type of planets beyond the snow line (Suzuki et al. 2016).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 −4 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(Suzuki et al. , 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.It is important to notice that the Suzuki et al. (2016) sample had only 30 exoplanets, which makes the expansion of this sample 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 (R.A., decl.)= (17 h 53 m 41 s .64,−29 • 48 27 .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.8m 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 nm 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 re-reduction 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 Mauna Kea in Hawaii, United States, also observed the event in the SDSS1 i-band filter.From 2020 March to July, CFHT provided 1-2 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 two-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 Korea Microlensing Telescope Network (KMTNet) group as KMT-2020-BLG-0579.The KMTNet collaboration monitored this event with the 1.6m 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 Interamerican 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 datasets 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.

LIGHT-CURVE MODELS
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 HJD = [9040.7,9041.5]2observed 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 propose the lack of evidence for a binary-source model -with single lens (1L2S) or binary lens (2L2S) -after investigating the light curve using the method described in Bennett et al. (2018).

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 χ 2 minima.

Fit Parameters
The parameters of our model are: the Einstein crossing time (t E ); the time at which the separation of lens and source reaches the minimum (t 0 ); the minimum angular separation between source and lens as seen by the observer (u 0 ); 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 0 and u 0 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 0 and s, are normalized by the angular Einstein radius θ E , defined by 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, 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 χ 2 .The total flux F i (t) for time t for instrument i can be written as: 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 0 , u 0 , 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 Earth's 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.

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 0 , u 0 , 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  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 −5 , 10 −1 ] 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 χ 2 , which indicates a model with q ∼ 10 −4 .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 χ 2 -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.2L1S), and, in the diagonal, the 1dimension Probability Density Function (PDF) of each parameter.The 68.3% (1σ), 95.5% (2σ), and 99.7% (3σ) confidence intervals are shown in dark blue, median blue, and light blue, respectively, in the posterior distribution plots.In the PDF plots, the black dot points out the median, and the thin line marks the 1σ confidence interval.For the separation s PDF plot, the additional red and yellow dots and lines also point out the median and the 1σ confidence interval, but now for each of the two regions, the close-separation (s < 1) in red and the wide-separation (s > 1) in yellow.
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.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 from our algorithm jumping back and forth in between both wide and close solutions.Even though there seems to exist two regions in each butterfly wing when considering the 1σ interval, the chi-square difference is not big enough to create a barrier that could separate them into four defined regions.Figure 3a and Figure 3b illustrate the close-separation (s < 1) and the wideseparation (s > 1) topology, respectively.We can see in Figure 3 that, although the topology of the caustics for each solution has a different shape, the source trajectory crosses parts with slightly similar shapes of both caustics.An (2021) argues that local degeneracies of lensing features associated with caustic crossings can still exist in the planetary events even in situations in which the overall caustic shape may not be degenerate at all.As a consequence, 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 planet presence appears to start at HJD ≈ 9040.85,finishing at HJD ≈ 9041.30, while for s > 1, it appears to start at HJD ≈ 9040.82,finishing at HJD ≈ 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 different time.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 close-wide degeneracy solutions.One might wonder about the proximity of the source and the upper cusp at HJD ≈ 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 HJD ≈ 9041.71.It results in no apparent anomaly in the best-fit model in Figure 1.

Binary Source Possibility
The investigation of the possibility of two sources was encouraged when still considering using the KMTNet data.However, the evidence of binary source largely disappeared when systematic of KMTNet data was discovered.
Following the suggestion of the KMTNet collaboration, we remove the KMTNet data from the analysis, as explained in Section 2, and we attempt to fit binary source models to our data.We use the method described in Bennett et al. (2018) to look for, first, solutions with binary-source single-lens (1L2S) and, then, with binary-source binary-lens (2L2S) models.For easier comparison, we refer to the Point-Source Point-Lens (PSPL) model as PSPL/1L1S, the best-fit model as 2L1S s > 1, and its degenerate model as 2L1S s < 1 (Section 3.1).
The Bennett et al. (2018) method allows us to add a second source with different brightness and color from the first source.We find the best-fit for the 1L2S model to have a fit χ 2 = 14862.24,with a lens-source proper motion of µ rel,G = 0.44 ± 0.05 mas yr −1 , which is unusually small.With ∆χ 2 1L2S−2L1S = 42.53 and a small relative proper motion, our results indicate a strong preference toward the 2L1S s > 1 model.Therefore, the 1L2S model is ruled out because it suggests that a model with a second source and only one lens reduces the quality of the fitting to the data.
For the binary-source binary-lens, 2L2S, model, we include the binary-lens parameters.We obtain a fit χ 2 = 14807.54,which is smaller than the one for the 2L1S s > 1 model (∆χ 2 2L1S−2L2S = −12.17).As expected, due to a larger number of parameters, a model with 2L2S may improve the χ 2 , so an extra criterion is necessary to evaluate the significance of this small improvement.Therefore, we also calculate the Bayesian Information Criterion (BIC) (Schwarz 1978) indexes for the two models.The difference in BIC index is ∆BIC 2L2S−2L1S = 35.85.The improvement of the quality of the fit is, then, proved to be insignificant when using the BIC index.Moreover, it should be noted that a binary source event requires special alignment of the sources, being generally unlikely to be detected.Therefore, we continue our analysis considering only the 2L1S model.The differences in χ 2 and BIC index for each model compared to our best-fit model is displayed in Table 2.

PHOTOMETRIC CALIBRATION AND SOURCE PROPERTIES
The source angular radius θ * is not explicitly obtained from the lightcurve models presented in Section 3, yet it can be empirically derived if we know the de-reddened 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 a standard calibration procedure described in Bond et al. (2017): where R MOA is the MOA-Red filter and V MOA is the MOA-V filter.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.
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 bulge3 , 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 −0.078 and dereddened color (V − I) S,0 = 0.788 +0.095 −0.096 (in Table 3) of our source star.
Finally, we determine the source size: 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 θ * = 1.15 +0.13 −0.12 µ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 , and ( 7) where the lens-source relative parallax π rel is defined as 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 θ E = 0.133 +0.019 −0.018 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: The lens-source proper motion is µ rel,G = 2.88 +0.42 −0.40 mas 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 K-band by using the color transformations presented in Kenyon & Hartmann (1995).The source magnitude is K S,0 = 16.45 +0.22  −0.25 without extinction.By adding the extinction A K = 0.2195 (Nishiyama et al. 2009;Gonzalez et al. 2012), we obtain K S = 16.67 +0.22   −0.25   as the magnitude of the source in K-band.
The calculated source and lens-source properties described in this section are in Table 3.

PHYSICAL PROPERTIES OF THE LENS
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  2014) Galactic model using a power-law stellar mass function with a prior uniform in M , dP ∝ dM , meaning all stars have an equal planet-hosting probability.The 1σ and 2σ confidence intervals (i.e., 68.3% and 95.5%) are represented by dark blue and median blue, respectively, with the median marked as a solid black line.For comparison, the source magnitude is shown as a green line in the lens magnitude graphs.
proportion to the host star mass, dP ∝ M dM .Johnson et al. (2007Johnson et al. ( , 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 planet = 11.3 +19.2 −6.9 M ⊕ for the planet mass, M host = 0.23 +0.39  −0.14 M for the host mass, D L = 7.9 +1.0 −1.0 kpc for the distance to the lens, a ⊥ = 1.11 +0.23  −0.20 AU for the project separation, a 3D = 1.35 +0.75 −0.32 AU for the deprojected separation, and I L = 26.0+2.4  −2.9 and K L = 22.2 +1.9 −2.3 for the lens magnitude.When assuming the mass function prior proportional to M , β = 1 ⇒ dP ∝ M dM , the planet mass estimation is m planet = 25 +22 −15 M ⊕ , the host mass is M host = 0.53 +0.42 −0.32 M , the distance to the lens is D L = 8.3 +0.9 −1.0 kpc, the projected separation is a ⊥ = 1.17 and K L = 20.4+2.0 −2.3 .We report the probability distributions considering both priors in Table 4.
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 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.

DISCUSSION
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.39  −0.31 × 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 K -band, aiming to anticipate results for high-angular resolution follow-up observations.The source magnitude with added extinction was computed to be K S = 16.67 +0.22  −0.25 .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 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.39  −0.31 × 10 −4 seems to be exactly where previous core accretion theories predicted a Neptune desert (Ida & Lin 2004;Mordasini et al. 2009).On the other hand, it also lies exactly in the peak of the planet-to-star mass ratio distribution measured by gravitational microlensing (Suzuki et al. 2018).Once more, the Suzuki et al. (2016) contained only 30 planets and, to better understand this contradiction, 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.
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.39  −0.31 × 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(Johnson et al. , 2010) ) and direct detection (Nielsen et al. 2019) surveys.However, both of these analyses considered planet hosting probabilities at fixed planet mass instead of fixed mass ratio.Since lower mass ratio planets are more common (Suzuki et al. 2016), for mass ratios ≥ 10 −4 , the demonstration that higher mass hosts are more likely for a fixed mass ratio is stronger than the same claim at fixed planet mass.If the planet hosting probability was independent of host mass at fixed mass ratio, the hosting probability would be higher for higher mass hosts at fixed planet mass, since this implies a larger mass ratio, q, for the lower mass hosts (as long as ≥ 10 −4 ).

CONCLUSION
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.39  −0.31 × 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 planet = 11.3 +19.2 −6.9 M ⊕ , a hoststar with mass of M host = 0.23 +0.39  −0.14 M and magnitude K L = 22.2 +1.9 −2.3 , located at a distance D L = 7.9 +1.0 −1.0 kpc.Under the assumption that the hosting-probability scales in proportion to the stellar mass, we have estimated m planet = 25 +22 −15 M ⊕ , M host = 0.53 +0.42  −0.32 M , K L = 20.4+2.0 −2.3 , and D L = 8.3 +0.9 −1.0 kpc.The previous results from RV and direct imaging (Johnson et al. 2007(Johnson et al. , 2010;;Nielsen et al. 2019) considered planet hosting probabilities for fixed planet mass, whereas the MOA-II exoplanet microlens statistical analysis considers 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.

Figure 1 .
Figure 1.The best-fit model for the MOA-2020-BLG-135 light curve in a time (HJD = HJD − 2, 450, 000) vs. magnification plot.The MOA-Red, MOA-V, and CFHT-i data are shown in brown, violet, and blue, respectively.The best single source planetary model 2L1S (i.e., two lenses and one source) s > 1 is our best-fit model (see Section 3) and is displayed as a black solid line, while the Point-Source Point-Lens (PSPL) model is displayed as a dashed gray line.The main panel shows all the photometric data sets overlaid with the best-fit model.The upper left side of the Figure, zooming near the peak of the event, shows the planetary perturbation, which is well covered by the MOA-Red data set, and the corresponding residuals.The lower panel shows the residuals of the best-fit model for each instrument -see colored legend.

Figure 2 .
Figure2.The marginalized posterior distributions for our MCMC runs, correlating the parameters for the binary lens single source model (2L1S), and, in the diagonal, the 1dimension Probability Density Function (PDF) of each parameter.The 68.3% (1σ), 95.5% (2σ), and 99.7% (3σ) confidence intervals are shown in dark blue, median blue, and light blue, respectively, in the posterior distribution plots.In the PDF plots, the black dot points out the median, and the thin line marks the 1σ confidence interval.For the separation s PDF plot, the additional red and yellow dots and lines also point out the median and the 1σ confidence interval, but now for each of the two regions, the close-separation (s < 1) in red and the wide-separation (s > 1) in yellow.
(a) Caustic geometry of the best-fit model for s < 1, due to a close-separation planet.(b) Caustic geometry of the best-fit model for s > 1, due to a wide-separation planet.

Figure 3 .
Figure 3.The caustic geometry and the source-lens trajectory.The caustic is represented in blue.The black straight solid line shows the source-lens trajectory and the arrow shows the direction of the source-lens relative motion.The source size is displayed as a black circle at its position at HJD = 9040.49(event peak), HJD = 9040.85for s < 1 and HJD = 9040.82for s > 1(source starts crossing the caustic), and HJD = 9041.00for s < 1 and HJD = 9041.34for s > 1 (source exiting the caustic).For (b), the source at HJD = 9041.71 is also displayed to show when it is closest to the upper cusp.

Figure 4 .
Figure 4. (V -I, I) Color-magnitude diagram of the stars in the OGLE-III catalog within 90 of MOA-2020-BLG-135.The black dots are the stars from the OGLE-III catalog, the blue dot indicates the source magnitude and color for the best-fit model (s > 1), and the red circle indicates the red clump giant centroid.For comparison, we added the green dots showing the Hubble Space Telescope CMD from Holtzman et al. (1998) shifted to the bulge distance and relevant extinction derived in Section 4. The source star is probably in a subgiant phase.

Figure 5 .
Figure 5. Lens system properties derived from the Bennett et al. (2014) Galactic model using a power-law stellar mass function with a prior uniform in M , dP ∝ dM , meaning all stars have an equal planet-hosting probability.The 1σ and 2σ confidence intervals (i.e., 68.3% and 95.5%) are represented by dark blue and median blue, respectively, with the median marked as a solid black line.For comparison, the source magnitude is shown as a green line in the lens magnitude graphs.

Figure 6 .
Figure 6.Lens system properties derived from the Bennett et al. (2014) Galactic model using a power-law stellar mass function under the assumption that the hosting probability scales in proportion to the stellar mass, dP ∝ M dM .The 1σ and 2σ confidence intervals (i.e., 68.3% and 95.5%) are represented by dark blue and median blue, respectively, with the median marked as a solid black line.For comparison, the source magnitude is shown as a green line in the lens magnitude graphs.

Table 1 .
Best-fit Model Parameters for MOA-2020-BLG-135, and Corresponding Averages from the Posterior Distribution version of the Metropolis algorithm, as implemented byBennett (

Table 2 .
Comparison Between Microlensing Models

Table 3 .
Source and lens-source properties