Kojima-1Lb Is a Mildly Cold Neptune around the Brightest Microlensing Host Star

We conducted photometric monitoring observations of Kojima-1 using 13 ground-based telescopes distributed around the world through the optical (g, r, i, z_s, B, V, R, and I) and near-infrared (K_s) bands, as listed in Table 1. The photometric follow-up campaign started on 2017 October 31 and lasted for 76 days until the source's brightness well returned to the original state. The number of observing nights, median observing cadence after removing outliers and time-binning, and median photometric error of each instrument are appended to Table 1. We note that we triggered the follow-up campaign without knowing the presence of the planetary anomaly, which was first reported on 2017 November 8 (Nucita et al. 2017). Also, we did not change any observing cadences after the report of the anomaly detection because (1) the anomaly had already finished at the time of the report and therefore no further follow-ups were required for the anomaly itself, and (2) from the beginning, we intended to follow up the event as much as possible until the end of the event, no matter whether a planetary anomaly was detected around the peak or not, to search for new planetary signals. On the other hand, we would have terminated our follow-up campaign by the end of 2017 if the planetary anomaly was not detected, and we extended the campaign for ∼2 weeks in reaction to the anomaly detection hoping to place a better constraint on the microlensing light-curve model. We will reflect this point in the calculation of the planet detection efficiency in Section 6.2. We further note that the data from CBABO and SL in the list were also used in Nucita et al. (2018); however, we rereduced them with our own photometric pipeline in order to investigate the possible systematics in these data (see below for CBABO and Section 3.3 for SL).

All of the data were corrected for bias and flat-field in a standard manner. To extract the light curves of the event, aperture photometry was performed using a custom pipeline (Fukui et al. 2011) for the data sets of MuSCAT, MuSCAT2, ISAS, OAOWFC, CBABO, COAST, SL, and MITSuME; IRAF/APPHOT⁴⁶ for Araki; SExtractor (Bertin & Arnouts 1996) for PROMPT-8; and AIJ (Collins et al. 2017) for OAR and WCO, and a differential image analysis using the ISIS package⁴⁷ (Alard & Lupton 1998; Alard 2000) was performed for the data set of DEMONEXT. In the case of aperture photometry, comparison stars are carefully selected for each data set depending on the field of view, so that systematics arising from intrinsic variabilities of the comparison stars are minimized.

On the raw images of CBABO obtained on 2017 October 31, the flux counts of the target star were close to the saturation of CCD and affected by the CCD nonlinearity. We corrected this effect by constructing a pixel-level nonlinearity-correction function using a seventh-order polynomial by minimizing the dispersion of the aperture-integrated light curve of a similar-brightness star in the same field of view (TYC 1849-1592-1).

The observed light curves are shown in Figure 1 in magnification scale. While we confirmed the planetary feature around the peak in the data sets of COAST, CBABO, and SL, we did not detect any additional anomaly in the light curves.

Zoom In Zoom Out Reset image size

A high-resolution spectrum was taken in the wavelength range of 4990–7350 Å using the NAOJ 188 cm telescope in Okayama, Japan, and the High Dispersion Echelle Spectrograph (HIDES; Kambe et al. 2013) on 2017 November 1.6 UT. Two exposures were obtained in the high-efficiency mode (HE mode; R ∼ 55,000) with exposure times of 23 and 20 minutes. The data reduction (bias subtraction, flat-fielding, spectrum extraction, and wavelength calibration) was performed using the IRAF echelle package in a standard manner. The signal-to-noise ratio (S/N) of the obtained spectrum is approximately 20–30.

Low-resolution spectra (R ∼ 500) were taken on 2017 November 3 and 2018 January 3 using the FLOYDS spectrograph mounted on the Las Cumbres Observatory (LCO) 2 m telescope on Haleakala, Hawaii.⁴⁸ The spectral range is about 3200–10000 Å. Each spectrum was taken with 1000 s exposure with the 12 slit. Both spectra were obtained in similar sky conditions, but due to the different magnification at the time of exposure (8.34 and 1.04), both images were obtained with different S/Ns, a range of [50, 250] and [20, 90], respectively. Both 1D spectra were extracted using the FLOYDS pipeline.⁴⁹

High-resolution images of the event object were obtained using the Keck telescope and NIRC2 instrument on 2018 February 5. Using the narrow camera (pixel scale of 9.94 mas pixel⁻¹), 10 dithered images were obtained in the K_s band with the NGS mode, each with an exposure time of 2 s and three coadds. The median FWHM of the adaptive optics (AO)–guided stellar point-spread function was 006. The raw images were median-combined after bias flat correction, sky subtraction, and stellar position alignment. The combined image and a 5σ contrast curve are shown in Figure 2. We found no contaminating sources brighter than K_s = 21 within the image.

Zoom In Zoom Out Reset image size

To derive the physical parameters of the lens system, we fit the light curves with a binary-lens microlensing model. The model calculates the magnification of the source star as a function of time, A(t), which is expressed by the following parameters: the time of the closest approach of the source to the lens centroid, t₀; the Einstein radius crossing time, t_E; the source-lens angular separation at time t₀ in units of the angular Einstein radius (θ_E), u₀; the mass ratio of the binary components, q; the sky-projected separation of the binary components in units of θ_E, s; the angle between the source trajectory and the binary-lens axis, α; the angular source radius in units of θ_E, ρ; and the microlens parallax vector, ${{\boldsymbol{\pi }}}_{{\rm{E}}}$. Here the direction of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ is the same as the direction of the source's proper motion relative to the lens, and the length of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$, ${\pi }_{{\rm{E}}}\equiv \sqrt{{\pi }_{{\rm{E}},{\rm{N}}}^{2}+{\pi }_{{\rm{E}},{\rm{E}}}^{2}}$, is equal to the ratio of 1 au to the projected Einstein radius onto the observer plane, where π_E,N and π_E,E are the north and east components of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$, respectively. The limb-darkening effect of the source star is modeled by the formula $I(\theta )=I(0)[1-{u}_{X}(1-\cos \theta )]$, where θ is the angle between the normal to the stellar surface and the line of sight, I (θ) is the stellar intensity as a function of θ, and u_X is a coefficient for filter X. The observed flux in the ith set of instrument and band at time t is expressed by the linear function ${F}_{i}(t)=A(t)\times {F}_{s,i}+{F}_{b,i}$, where F_s,i and F_b,i are the unmagnified source flux and blending flux, respectively, in the ith data set. Note that the effect of the orbital motion of the planet is not considered in the final analysis because it was not significant in the first trials.

The initially estimated uncertainties of individual data points are rescaled using the formula

$$\begin{eqnarray}&&{\sigma }_{i}^{\prime} =k\sqrt{{\sigma }_{i}^{2}+{e}_{\min }^{2}},\end{eqnarray} \tag{ 4 }$$

where σ_i is the initial uncertainty of the ith data point in magnitude, and k and e_min are coefficients for each data set. Here the term e_min represents systematic errors that dominate when the flux is significantly increased. The k and e_min values are adjusted so that the cumulative χ² distribution for the best-fit binary-lens model including the parallax effect sorted by magnitude is close to linear and ${\chi }_{\mathrm{red}}^{2}$ becomes unity. This process is iterated several times.

In addition, we quadratically add 0.5% in flux to each flux error for the data points that lie within the anomaly, taking into account the possible intrinsic variability of the target and/or comparison stars. This additional error is important to properly estimate the uncertainties of the model parameters, in particular s, ρ, and π_E, which we find are sensitive to this anomaly part and can be biased by even a small systematics of the level of 0.5% in flux.

To save computational time, we restrict the data sets for a light-curve fitting to the ones with relatively high photometric precision with sufficient time coverage and/or unique coverage in time or wavelength, specifically, the data sets of MuSCAT, MuSCAT2, Araki, ISAS, OAOWFC, CBABO, and COAST. To supplement our data, we also use the V-band light curve from the All-Sky Automatic Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017; data are extracted from their website⁵⁰ for the period of 7967 < HJD–2,450,000 < 8123), which covered the entire event with an average cadence of several per night, and the V-band light curve capturing the declining part of the anomaly obtained at the R. P. Feynman Observatory (FO) by Nucita et al. (2018).

We note that although the SL data set includes the earliest data points among all of the follow-up observations partly overlapping with the FO data set (HJD–2,450,000 ∼ 8058.5), we have not included it in our light-curve modeling for the following reasons. First, when we fit the light curves including this data set, we found that the data points of this data set in the anomaly part have a small systematic trend against the best-fit model. Second, we also found that the F_s and F_b values from this data set, calibrated to standard photometric systems, were discrepant with those from the other same-band data sets at the 2σ level,⁵¹ even using only the data points that overlap with the FO. Because light-curve models are sensitive to the data points in the anomaly part, even a 2σ level systematics could cause tension in the derived parameters.

The light curves are fitted with a binary-microlensing model using a custom code that has been developed for the Microlensing Observations in Astrophysics (MOA) project (Sumi et al. 2010), in which the posterior probability distributions of the parameters are calculated by the Markov chain Monte Carlo (MCMC) method. Note that the light curves are also independently analyzed using the pipeline PyLIMA (Bachelet et al. 2017), a code developed by Bennett (2010), and the modeling platform RTModel⁵² (Bozza et al. 2018) for sanity check.

We first fit the light curves with a binary-lens model without the microlens parallax effect (static model), fixing ${\pi }_{\mathrm{EE}}$ and π_EN at zero, to compare with the result of Nucita et al. 2018, in which this effect was not taken into account. The median value and 1σ confidence interval of the posterior probability distributions of the parameters are listed in Table 2. We recover the two degenerate models found by Nucita et al. (2018; models a and b), in which only s is slightly different and all the other parameters are almost identical between the two models. The best-fit χ² values are almost the same between the two models, namely, 2557.5 and 2557.4 for models a and b, respectively, for the degrees of freedom (dof) of 2578. In Table 2, we report the values derived only for model b for all parameters except for s, and hereafter, we will discuss them along with this model unless otherwise described.

Table 2.  Best-fit Parameter Values of Binary-lens Microlensing Models

Parametera	Unit	Nucita et al. (2018)	Static	Parallax	Parallax	Parallax
					w/θE Prior	w/θE, Φπ Priors
t0	HJD	0.75 ± 0.01	0.7353 ± 0.0076	0.7395 ± 0.0073	0.7396 ± 0.0073	0.7403 ± 0.0074
	−2,458,058
tE	days	26.4 ± 0.9	27.44 ± 0.07	27.19 ± 0.15	27.18 ± 0.14	27.25 ± 0.09
u0	10−2	9.3 ± 0.1b	${8.858}_{-0.034}^{+0.031}$	8.925 ± 0.043	8.927 ± 0.042	8.935 ± 0.038
q	10−4	1.1 ± 0.1	${1.058}_{-0.074}^{+0.068}$	${1.075}_{-0.073}^{+0.066}$	${1.031}_{-0.084}^{+0.078}$	${1.027}_{-0.084}^{+0.078}$
s (model a)		0.935 ± 0.004	${0.9207}_{-0.0040}^{+0.0045}$	${0.9204}_{-0.0038}^{+0.0040}$	$0.9263\,\pm \,0.0018$	$0.9264\,\pm \,0.0018$
s (model b)		0.975 ± 0.004	${0.9944}_{-0.0046}^{+0.0041}$	${0.9941}_{-0.0045}^{+0.0042}$	$0.9874\,\pm \,0.0018$	$0.9873\,\pm \,0.0018$
α	rad	$4.767\pm 0.007$ b	$4.7594\,\pm \,0.0030$	$4.7610\,\pm \,0.0030$	$4.7604\,\pm \,0.0028$	4.7604 ± 0.0028
ρ	10−3	6.0 ± 0.8	${3.2}_{-1.3}^{+0.9}$	${3.2}_{-1.3}^{+0.9}$	$4.568\,\pm \,0.070$	4.567 ± 0.071
${\pi }_{{\rm{E}},{\rm{E}}}$		⋯	⋯	${0.071}_{-0.064}^{+0.072}$	${0.0693}_{-0.063}^{+0.070}$	${0.143}_{-0.053}^{+0.061}$
${\pi }_{{\rm{E}},{\rm{N}}}$		⋯	⋯	$0.17\,\pm \,0.45$	$0.19\,\pm \,0.45$	$-{0.33}_{-0.14}^{+0.12}$
${\chi }_{\min }^{2}$ /dof		⋯	2557.4/2578	2543.0/2576	2546.5/2577	2550.7/2578
${\pi }_{{\rm{E}}}$		⋯	⋯	${0.34}_{-0.20}^{+0.34}$ c	${0.35}_{-0.20}^{+0.34}$ c	${0.36}_{-0.13}^{+0.16}$
${\theta }_{* }$	μas	⋯	8.59 ± 0.06	8.65 ± 0.06	8.63 ± 0.06	8.63 ± 0.06
${\theta }_{{\rm{E}}}$	mas	1.45 ± 0.25	${2.63}_{-0.58}^{+1.77}$	${2.68}_{-0.59}^{+1.87}$	1.890 ± 0.032	1.890 ± 0.032

Notes.

^aThe values for the two models (a and b) are basically identical except for s, for which both values are presented. Only the values for model b are presented for the other parameters. ^bFor ease of comparison, we multiply the u₀ and increment α reported in the literature by −1 and π, respectively. The geometry is identical to this transformation. ^cBecause ${\pi }_{{\rm{E}},{\rm{E}}}$ and π_E,N take both positive and negative values, the median value of π_E does not coincide with $\sqrt{{\left\langle {\pi }_{{\rm{E}},{\rm{E}}}\right\rangle }^{2}+{\left\langle {\pi }_{{\rm{E}},{\rm{N}}}\right\rangle }^{2}}$, where $\left\langle {\pi }_{{\rm{E}},{\rm{E}}}\right\rangle$ and $\left\langle {\pi }_{{\rm{E}},{\rm{N}}}\right\rangle$ are the median values of π_E,E and π_E,N, respectively.

Download table as:  ASCIITypeset image

Our derived values are consistent with those of Nucita et al. (2018) within 2σ for all parameters except for u₀, s, and ρ, for which the discrepancy can be attributed to the following differences between our and their data sets: (1) we correct the detector's nonlinearity effect in the CBABO data set, (2) we omit the SL data set from our modeling due to apparent systematics, and (3) we have a larger number of data points with a longer baseline.

3.5.1. Without Informative Prior

To search for a signal of the parallax effect, we fit the light curves letting π_EE and π_EN be free, first without any informative priors. The derived values and uncertainties are reported in Table 2. From this fit, we marginally detect a nonzero π_E value of ${0.34}_{-0.20}^{+0.34}$. However, the χ² improvement of the best-fit parallax model over the static model is 14.4, which is not significant enough to claim a detection of the parallax signal, given that the Bayesian information criterion (BIC $\equiv \,{\chi }^{2}+k\mathrm{ln}{N}_{\mathrm{data}}$, where k is the number of free parameters and N_data = 2615 is the number of data points) for the parallax model is larger (worse) than the static model by 1.3.

We also check where the marginal parallax signal comes from. In the top panel of Figure 3, we show the magnitude differences between the best-fit static and parallax models for individual data sets, which indicate that the largest difference arises around ∼20 days before the peak, yet the difference is at most at the ∼10 mmag level. On the other hand, in the bottom panel of Figure 3, we show the difference of cumulative χ² between the two models as a function of time. This plot indicates that most of the χ² improvements come from only two epochs of the MuSCAT data (from three different bands), where the model magnitudes differ by only ∼1 mmag. Thus, the likely origin of the parallax signal is due to systematics in the data at these two epochs, which might arise from the instrument, variability of atmospheric transparency, and/or stellar activity. Therefore, the observed marginal signal of the parallax effect should be treated with caution. Nevertheless, the data still allow us to place an upper limit on π_E (Section 3.5.3) and constrain the direction of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ (Section 3.5.4).

Zoom In Zoom Out Reset image size

The result that a significant parallax signal is absent is consistent with the result of Dong et al. (2019), who also did not detect a significant parallax signal from a single-lens model fit (for the "luminous-lens" case in their paper). Dong et al. (2019) described the reasons why the parallax signal in this event is not obvious, which are summarized as follows: (1) the event is quite short compared to a year, (2) it lies quite close to the ecliptic plane, (3) it peaked only 5 weeks⁵³ before opposition, and (4) the lens-source relative proper motion points roughly south. The combination of these factors weakens the parallax signal in the light curve by a factor of ∼10 compared to the most favorable case (Dong et al. 2019).

3.5.2. With Informative Prior on θ_E

3.5.3. Upper Limit on π_E

From the VLTI observation, Dong et al. (2019) also constrained the direction of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ (Φ_π) into two directions, 1935 ± 04 and 1567 ± 04 from north to east (for the luminous-lens model). To put an upper limit on π_E utilizing the prior information of Φ_π, we draw χ² maps on a grid of π_E,E and π_E,N. We grid π_E,E and π_E,N by a grid size of 0.1 in the ranges of $-0.7\,\leqq \,{\pi }_{{\rm{E}},{\rm{E}}}\lt 0.7$ and $-1.5\,\leqq \,{\pi }_{{\rm{E}},{\rm{N}}}\lt 1.5$ and fit the light curves using the θ_E prior while fixing π_E,E and π_E,N at each grid-point value. In the left panel of Figure 4, we show Δχ² maps on the π_E,E–π_E,N plane calculated from all data sets, where Δχ² is the difference of χ² between each grid point and (π_E,E, π_E,N) = (0, 0). The minimum-χ² (darkest red) region is not coincident with the two solutions of Φ_π (indicated by cyan lines), probably due to the systematics in the light curves discussed before. Note that the reason why the negative Δχ² region is elongated almost along the π_E,N direction (only π_E,E is well constrained) is that the direction of Earth's acceleration is almost parallel to the direction of π_E,E.⁵⁴ On the other hand, the right panel of the same figure shows a Δχ² map that is calculated only using the χ² values from the ASAS-SN data set, which covers the region where the parallax signal is maximized and is thus robust for a parallax signal against the systematics. In this map, although the minimum-χ² region is not localized, the intersection between the Φ_π solutions and some Δχ² contour can still be used to put an upper limit on π_E. The contour of Δχ² = 9 (white) intersects with the Φ_π ∼ 1567 and 1935 lines (cyan) at the grid points that correspond to π_E = 1.1 and 0.5, respectively. We conservatively adopt 1.1 as a 3σ upper limit on π_E.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.5.4. On the Direction of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$

As will be discussed in Section 5.2.2, under the condition of π_E < 1.1, it is most likely that the blending flux detected in the light curves comes from the lens star independently on the Φ_π value, and this lens flux allows us to derive the mass of the lens star to be ${M}_{L}={0.590}_{-0.051}^{+0.042}{M}_{\odot }$. This lens mass, combined with θ_E, predicts the π_E value using the relation

$$\begin{eqnarray}&&{\pi }_{{\rm{E}}}=\displaystyle \frac{{\theta }_{{\rm{E}}}}{\kappa {M}_{L}},\end{eqnarray} \tag{ 5 }$$

where $\kappa \equiv 4G/{c}^{2}$, G is the gravitational constant, and c is the speed of light. This gives ${\pi }_{{\rm{E}}}={0.39}_{-0.03}^{+0.04}$, which is indicated by magenta solid (median) and dotted (1σ boundary) contours in Figure 4. In the Δχ² map for all data sets (left panel of Figure 4), the Δχ² value at the grid point that satisfies both π_E ∼ 0.39 and Φ_π ∼ 1567 is −16, which is smaller than the counterpart that satisfies both π_E ∼ 0.39 and Φ_π ∼ 1935 by 40. This χ² difference nominally rules out the Φ_π = 1935 solution.

This outcome, however, could be affected by systematics in the light curves. To test this possibility, we also check the Δχ² map calculated using only the χ² values from the ASAS-SN data set (right panel of Figure 4). We find that the Φ_π = 1567 solution is preferred over the other solution with a χ² improvement of ∼5, which, although marginal, supports the outcome obtained from all data sets.

Considering the above evidence, we adopt the ${{\rm{\Phi }}}_{\pi }=156\buildrel{\circ}\over{.} 7$ solution for further analysis. To derive the final posteriors of the parameters, taking into account correlations between the parallax parameters (π_E,E and π_E,N) and others and using all of the informative prior information, we rerun the MCMC analysis, letting π_E,E and π_E,N be free and imposing priors on θ_E and Φ_π with Gaussian distributions of θ_E = 1.883 ± 0.014 mas and Φ_π = 1567 ± 04. The results are reported in Table 2, and the caustics for the models a and b, along with the source trajectory and northeast direction, are shown in Figure 5. We note that if the other solution of Φ_π is adopted, then the light-curve fit gives slightly larger values of the blending flux, leading to an ∼10% increase of M_L. This, however, does not change the conclusion of this paper much. This Φ_π value can be confirmed in the future by directly measuring the lens-source relative position from high spatial resolution images.

The spectroscopic properties of the source star are initially estimated from the HIDES spectrum in the wavelength region of 5000–5900 Å. Note that the spectrum in longer wavelengths is not used to avoid a significant fringe effect. Because the spectrum was taken at a time when the source was magnified by a factor of 10, the flux contamination from other objects into the source's spectrum is negligibly small, with a fraction of less than 0.4% in this wavelength range. We also note that the spectrum does not show any sign of a companion star, i.e., a split of lines due to differential radial velocity. Using the spectral fitting tool SpecMatch-Emp (Yee et al. 2017), which matches an observed spectrum with empirical spectral libraries, we estimate the stellar effective temperature, radius, and metallicity to be ${T}_{\mathrm{eff}}=6303\pm 110$ K, R_S = 1.56 ± 0.25 R_⊙, and [Fe/H] = −0.11 ± 0.08, respectively. This result indicates that the source star is a main-sequence late-F dwarf.

The two LCO spectra were taken at magnifications of A₁ = 8.34 and A₂ = 1.04, with which the flux contamination from the lens star, in particular for the wavelength of ≳700 nm, is not negligible. Nevertheless, we can extract the source spectrum from the observed spectra using the equation ${f}_{s,\lambda }=({f}_{1,\lambda }-{f}_{2,\lambda })/({A}_{1}-{A}_{2})$, where f_1,λ and f_2,λ are the fluxes at the wavelength λ in the first- and second-epoch spectra, respectively. We correct the interstellar extinction in the source spectrum and compare it with empirical spectral templates of Kesseli et al. (2017), as shown in Figure 6, finding that the source's spectral type is F5V ± 1 subtype. This result is consistent with that obtained from the HIDES spectrum.

Zoom In Zoom Out Reset image size

The interstellar extinction toward the source star is initially estimated using the Gaia DR2 (Gaia Collaboration et al. 2016, 2018), in which the trigonometric parallax (π) and extinction in the Gaia band (A_G) are both recorded for a subset of relatively bright and nearby stars. Although the uncertainties of individual A_G values are large, an ensemble of A_G can be used to estimate the averaged A_G value in the field because the uncertainties are dominated by statistical errors (Gaia Collaboration et al. 2018).

First, from the Gaia DR2, we extract stars that lie within 30' of the source position, have records of both π and A_G, and have π > 0.5 mas with a fractional uncertainty of less than 20%. Next, all of the data are divided by distance into bins with a width of 50 pc. The mean and 1σ error (standard deviation divided by the square root of the number of data points) for each bin are calculated, where the median 1σ error is ∼0.10. The binned data are then fitted with a fourth-order polynomial function of the distance, which gives

$$\begin{eqnarray}\begin{array}{rcl}{A}_{G} & = & -7.4918\times {10}^{-2}+3.6988\times {10}^{-3}D\\ & & -5.1142\times {10}^{-6}{D}^{2}+3.0569\times {10}^{-9}{D}^{3}\\ & & -6.4472\times {10}^{-13}{D}^{4},\end{array}\end{eqnarray} \tag{ 6 }$$

where D is the distance from the Earth. We plot the individual and binned A_G data along with the derived function in Figure 7. We also calculate the ratio of A_G to A_V, which is the extinction in the V band, to be 1.13, assuming the extinction law of Cardelli et al. (1989) with ${R}_{V}\equiv {A}_{V}/E(B-V)=3.1$.

Zoom In Zoom Out Reset image size

Although the trigonometric parallax of an object at the same coordinates as Kojima-1 was measured by Gaia to be 1.45 ± 0.03 mas, this value does not represent the true trigonometric parallax of the source star but is biased by the foreground lens star. Based on the multiband measurements of F_s and F_b, we estimate that the flux ratio of the lens to the source stars in the Gaia band is ∼5%, assuming that F_b comes entirely from the lens star (see Section 5.2.1). On the other hand, the Gaia DR2 data were acquired during the period between 3.3 and 1.4 yr before the peak of the event, which translates to lens-source separations of ∼83  and ∼35 mas, respectively. Because the image resolution of Gaia is 250 mas × 85 mas, this lens flux fully contaminated to the Gaia images, substantially changing its position relative to the source star. Therefore, it is not possible to estimate the effect of the lens-flux contamination on the measured parallax without knowing the respective times of the time series of Gaia astrometric data.

We instead estimate the distance (D_S) and angular radius (θ_*) of the source star using the spectral energy distribution (SED) as follows. First, we calibrate the source fluxes, F_s, in the g, r, i, and z_s bands of MuSCAT and MuSCAT2 to the SDSS g', r', i', and z' magnitudes, respectively. We also convert the F_s in the V band of ASAS-SN to the Johnson V magnitude and calibrate the F_s in the K_s band of OAOWFC to the 2MASS K_s magnitude (Table 3). The calibrated magnitudes are then converted into flux densities to create the SED. Next, we fit the SED with the synthetic spectra of BT-Settl (Allard et al. 2012) using the following parameters: the stellar effective temperature T_eff, radius R_S, metallicity [M/H], A_V to the source star A_V,S, and D_S. For a given set of R_S and [M/H], log surface gravity (log g) is calculated using an empirical relation of Torres et al. (2010), and from a set of T_eff, [M/H], and log g, a synthetic spectrum is created by linearly interpolating the grid models. The synthetic spectrum is then scaled by ${({R}_{S}/{D}_{S})}^{2}$ and reddened using a given A_V,S value and R_V = 3.1 to fit the observed SED. We perform MCMC to calculate the posterior probability distribution of each parameter using the emcee code (Foreman-Mackey et al. 2013). In the MCMC sampling, Gaussian priors are applied to the parameters T_eff, R_S, [M/H], and A_V,S by adding penalties to the χ² value as

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \sum _{\lambda }\displaystyle \frac{{\left({f}_{\mathrm{obs},\lambda }-{f}_{\mathrm{model},\lambda }\right)}^{2}}{{\sigma }_{{f}_{\mathrm{obs},\lambda }}^{2}}\\ & & +\sum _{i}\displaystyle \frac{{\left({X}_{i}-{X}_{i,\mathrm{prior}}\right)}^{2}}{{\sigma }_{{X}_{i,\mathrm{prior}}}^{2}},\end{array}\end{eqnarray} \tag{ 7 }$$

where ${f}_{\mathrm{obs},\lambda }$, ${\sigma }_{{f}_{\mathrm{obs},\lambda }}$, and ${f}_{\mathrm{model},\lambda }$ are the observed flux density, its 1σ uncertainty, and the model flux density, respectively, for a band λ, and X_i denotes one of the parameters among T_eff, R_S, [M/H], and A_V. For the priors of T_eff, R_S, and [M/H], the values derived from the HIDES spectrum are used, where [M/H] and [Fe/H] are assumed to be identical. As for A_V,S, the prior value is evaluated using Equation (6) for a given D_S, and 0.10 is taken as the 1σ uncertainty.

The derived median value and 1σ uncertainties of the parameters are reported in Table 3, and the posterior distributions are plotted in Figure 8. We derive the distance and angular radius of the source star to be D_S = 800 ± 130 pc and θ_* = 8.63 ± 0.06 μas, respectively, which are well consistent with the previous estimations of D_S = 700–800 pc (Nucita et al. 2018) and θ_* = 9 ± 0.9 μas (Dong et al. 2019).

Zoom In Zoom Out Reset image size

If θ_E, π_E, and D_S are all measured, one can solve for the total mass, M_L, and distance, D_L, of the lens system using the following formulae:

$$\begin{eqnarray}&&{M}_{L}=\displaystyle \frac{{\theta }_{{\rm{E}}}}{\kappa {\pi }_{{\rm{E}}}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{D}_{L}=\displaystyle \frac{{AU}}{{\pi }_{{\rm{E}}}{\theta }_{{\rm{E}}}+{\pi }_{S}},\end{eqnarray} \tag{ 9 }$$

where ${\pi }_{S}\equiv {AU}/{D}_{S}$. The masses of the host star and planet of the lens system are then calculated as ${M}_{L1}=1/(1+q){M}_{L}$ and ${M}_{L2}=q/(1+q){M}_{L}$, respectively, and the projected separation between the two lens components is derived by ${a}_{\mathrm{proj}}=s{\theta }_{{\rm{E}}}{D}_{L}$. The median and 1σ uncertainties of these parameters derived from the light-curve analysis using the θ_E and Φ_π priors (Section 3.5.4) are reported in Table 4, and the 68% and 95% confidence intervals of M_L1 and D_L are shown by blue dotted lines in Figure 9.

Zoom In Zoom Out Reset image size

Table 4.  Physical Parameters of the Lens System

Parameter	Unit	Nucita et al. (2018)	πE and θE and DS	Lens Flux	Lens Flux and θE and DS
Distance, DL	pc	∼380	${511}_{-80}^{+101}$	507 ± 74	$505\pm 47$
Stellar mass, ML1	M⊙	0.25 ± 0.18	${0.64}_{-0.19}^{+0.38}$	${0.590}_{-0.051}^{+0.042}$	0.586 ± 0.033
Stellar radius, RL1	R⊙	⋯	⋯	${0.599}_{-0.061}^{+0.056}$	⋯
Extinction, AV,L		⋯	⋯	0.95 ± 0.11	⋯
Metallicity, [Fe/H]	dex	⋯	⋯	−0.05 ± 0.20	⋯
Absolute Ks magnitude, ${M}_{{K}_{s}}$	mag	⋯	⋯	${5.05}_{-0.28}^{+0.33}$	⋯
Planetary mass, ML2	M⊕	9.2 ± 6.6	${21.8}_{-6.5}^{+12.9}$	20.0 ± 2.3	20.0 ± 2.0
Projected separation, aproj (model a)	au	∼0.5	${0.89}_{-0.14}^{+0.18}$	0.89 ± 0.13	0.88 ± 0.08
Projected separation, aproj (model b)	au	∼0.5	${0.95}_{-0.15}^{+0.19}$	0.95 ± 0.14	0.94 ± 0.09
Semimajor axis, acirca	au	⋯	${1.12}_{-0.25}^{+0.66}$	${1.10}_{-0.22}^{+0.63}$	${1.08}_{-0.18}^{+0.62}$

Note.

^aCalculated by merging the posteriors of models a and b.

Download table as:  ASCIITypeset image

However, as discussed in Section 3.5.1, the detection of π_E is marginal, and the signal is as weak as the level of systematics. Therefore, it is conservative not to rely on the π_E measurement to derive the lens parameters. In this case, we cannot uniquely solve for M_L1 and D_L but can only draw a relation between them, as shown by the gray shaded region in Figure 9.

5.2.1. Probabilities of Flux Contamination

From the light-curve fitting, we clearly detect the blending flux in the photometric aperture, F_b, in the optical and near-infrared bands from g through K_s. The F_b values in the g, r, i, z_s, V, and K_s bands are converted to the SDSS g', r', i', and z'; Johnson V; and 2MASS K_s magnitudes, respectively, as listed in Table 5.

Table 5.  Calibrated Magnitudes of the Blending Flux

Band	Magnitude
$g^{\prime}$	19.088 ± 0.337
V	17.760 ± 0.110
$r^{\prime}$	17.305 ± 0.122
$i^{\prime}$	16.382 ± 0.068
$z^{\prime}$	15.872 ± 0.051
Ks	13.728 ± 0.027

Download table as:  ASCIITypeset image

Generally, there are four possible sources that could contribute to the blending flux: the lens host, unrelated ambient stars, a companion to the source star, and a companion to the lens star. In the case of this event, however, the contribution from the ambient stars is negligible because the Keck AO image shows no stars with K_s < 21 mag in the sky area of 8'' × 8'' other than the target.

Following the method developed by Koshimoto et al. (2017) and N. Koshimoto et al. (2019 in preparation), we calculate the probabilities of all possible combinations of the other three sources that explain the observed blending flux, the Keck contrast curve (Figure 2), and the fact that the light curve shows no significant signal of a companion. In the calculation, we use the observed source and blending fluxes in the V, I, and K_s bands, where the fluxes in the I band are converted from those of i'- and z'-band magnitudes. We do not include stellar remnants. Using the posterior distribution from the MCMC calculation with the θ_E prior and the upper limit on π_E (<1.1), we calculate the probability distributions of the fraction of the lens flux to the total blending flux, f_L ≡ F_L/F_b, where F_L is the flux from the lens star. We find that the probability of f_L > 0.90 is 91.8%, which indicates that most of the blending flux most likely comes from the lens star. In the rest of the paper, we simply assume that the blending flux comes solely from the lens star. We note that the mass and distance of the lens star derived from the blending flux under the above assumption are well consistent with the constraint from θ_E and D_S (Section 5.1), supporting this assumption. There is still a small probability (8.2%) that more than 10% of the blending flux comes from a companion to the lens or source stars, which can be tested by direct imaging or spectroscopy of the lens star in the future.

5.2.2. Estimation of the Mass and Distance

With the assumption that the blending flux comes solely from the lens star, we can estimate the mass and distance of the lens star using the multiband blending flux. From an initial investigation, we find that the observed magnitudes and colors of the lens star are consistent with a main-sequence low-mass star. In estimation of the mass of low-mass stars, it is generally more reliable to use an empirical way rather than theoretical models (e.g., Boyajian et al. 2012). Therefore, to estimate a more accurate mass of the lens star, we adopt a mass–luminosity relation of Mann et al. (2019), which is a fully empirical and precise (2%–3% error on mass) mass–absolute-K_s relation for stars with a mass between 0.075 and 0.7 M_⊙, derived based on the apparent K_s magnitudes, trigonometric parallaxes, and dynamically determined masses of visual binaries. However, Mann et al. (2019) provided the relation only in the K_s band, with which alone the mass and distance of the lens star are degenerate for a given apparent K_s-band magnitude.

We therefore first solve for the distance and absolute K_s magnitude, ${M}_{{K}_{s}}$, from the apparent g'-, r'-, V-, i'-, z'-, and K_s-band magnitudes of the host star using empirical radius–metallicity–luminosity relations from Mann et al. (2015). They provided the relations based on spectroscopically measured effective temperatures, bolometric fluxes, metallicities, and trigonometric parallaxes of nearby M–K dwarfs in the form of

$$\begin{eqnarray}&&{R}_{* }=\sum _{i}^{n}{a}_{i}{M}_{\lambda }^{i}\times (1+f[\mathrm{Fe}/{\rm{H}}]),\end{eqnarray} \tag{ 10 }$$

where R_* is the stellar radius, M_λ is the absolute magnitude in the λ band, and a_i and f are coefficients. Because only the coefficients for the K_s band are provided in their paper, while they also collected apparent magnitudes in other bands, including the g', r', V, i', and z' bands, we derive the coefficients for these additional bands from the data sets of Mann et al. (2015) in the same way as they did for the K_s band (see the Appendix). We fit the observed magnitudes of the host star with a prediction calculated by

$$\begin{eqnarray}&&{m}_{\lambda ,\mathrm{calc}}={M}_{\lambda }+5{\mathrm{log}}_{10}({D}_{L}/10\mathrm{pc})+{A}_{\lambda ,L},\end{eqnarray} \tag{ 11 }$$

where λ is a given band, D_L is the distance to the lens in pc, and ${A}_{\lambda ,L}\equiv {A}_{V,L}\times {A}_{\lambda }/{A}_{V}$ is the extinction to the lens in the λ band. Note that M_λ is tied with the radius, R_L1, and metallicity, [Fe/H], of the lens star via Equation (10). Here we adopt A_λ/A_V = (1.223, 1.011, 0.880, 0.676, 0.485, 0.117) for λ = (g', r', V, i', z', K_s), calculated assuming R_V = 3.1.

We perform MCMC to derive the posterior distributions of D_L, R_L1, [Fe/H], and A_V,L using the emcee code (Foreman-Mackey et al. 2013). In this calculation, we evaluate the following χ² value:

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \sum _{\lambda =\{g^{\prime} ,r^{\prime} ,V,i^{\prime} ,z^{\prime} ,{K}_{s}\}}\displaystyle \frac{{\left({m}_{\lambda ,\mathrm{obs}}-{m}_{\lambda ,\mathrm{calc}}\right)}^{2}}{{\sigma }_{{m}_{\lambda ,\mathrm{obs}}}^{2}}\\ & & +\displaystyle \frac{{\left([\mathrm{Fe}/{\rm{H}}]-{\left[\mathrm{Fe}/{\rm{H}}\right]}_{\mathrm{prior}}\right)}^{2}}{{\sigma }_{{\left[\mathrm{Fe}/{\rm{H}}\right]}_{\mathrm{prior}}}^{2}}\\ & & +\displaystyle \frac{{\left({A}_{V,L}-{A}_{V,L,\mathrm{prior}}\right)}^{2}}{{\sigma }_{{A}_{V,L,\mathrm{prior}}}^{2}},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${m}_{\lambda ,\mathrm{obs}}$ and ${\sigma }_{{m}_{\lambda ,\mathrm{obs}}}$ are the observed magnitude and its 1σ uncertainty in the λ band, respectively; ${[\mathrm{Fe}/{\rm{H}}]}_{\mathrm{prior}}$ is a prior for [Fe/H]; and ${A}_{V,L,\mathrm{prior}}$ is a prior for A_V,L. Because our data alone do not put any meaningful constraint on [Fe/H], we impose a Gaussian prior with [Fe/H] = −0.05 ± 0.20, which is from the metallicity distribution of a nearby M dwarf sample (Gaidos & Mann 2014). We also take advantage of the extinction measurements of Gaia by applying Equation (6) to ${A}_{V,L,\mathrm{prior}}$ and 0.10 to ${\sigma }_{{A}_{V,L,\mathrm{prior}}}$ in the same way as for A_V,S. The derived posterior distributions of R_L1 and [Fe/H] are used to calculate the probability distribution of ${M}_{{K}_{s}}$ via Equation (10), which then gives the probability distribution of M_L1 via the mass–luminosity relation of Mann et al. (2019; Equation (2) of their paper where n = 5 is applied).

The derived median value and 1σ uncertainties of the parameters are presented in Table 4, and the posterior distributions of the parameters are plotted in Figure 10. The posterior distribution between D_L and M_L1 is also plotted in red in Figure 9. The derived D_L and M_L1 are well consistent with the constraints from the microlensing model (blue dotted contours and gray shaded region in Figure 9), while M_L1 is much better constrained by the lens flux.

Zoom In Zoom Out Reset image size

We derive the final values of M_L1 and D_L by combining the two posterior distributions, one from the microlens model (Section 5.1) and the other from the lens brightness (Section 5.2.2). For the microlens model, we use the posterior distribution of the M_L1–D_L relation derived from ${\theta }_{{\rm{E}}}$ and D_S instead of the posterior distribution of the M_L1 and D_L solution from ${\pi }_{{\rm{E}}}$, ${\theta }_{{\rm{E}}}$, and D_S, because the latter relies on the posterior distribution of ${\pi }_{{\rm{E}}}$, which could be affected by systematics (Section 3.5). Note that the posterior distribution from the lens flux and that from the microlens model can, in principle, be correlated because the blending flux that the former solution relies on was also derived using the microlens model. However, this effect is so small that these two distributions can be considered to be independent.

The combined posterior distribution is shown in green in Figure 9. As a result, we find that ${D}_{L}=505\pm 47$ pc and ${M}_{L1}=0.586\pm 0.033$ ${M}_{\odot }$; thus, the host star is a late-K/early-M boundary dwarf. The planetary mass is ${M}_{L2}\equiv {{qM}}_{L1}=20.0\,\pm 2.0$ ${M}_{\oplus }$, which is similar to the mass of Neptune (17.2 ${M}_{\oplus }$). The sky-projected separation between the planet and the host star is ${a}_{\mathrm{proj}}\equiv s{\theta }_{{\rm{E}}}{D}_{L}=0.88\,\pm \,0.08\,\mathrm{au}$ (model a) and $0.94\,\pm 0.09\,\mathrm{au}$ (model b), which is converted to the semimajor axis of ${a}_{\mathrm{circ}}={1.08}_{-0.18}^{+0.62}$ au, where a circular orbit and random orientation are assumed and the solutions of two models (model a and b) are merged.

Figure 11(a) shows the location of Kojima-1Lb in the plane between the mass and semimajor axis, along with the known exoplanets hosted by stars with masses similar to that of Kojima-1L (0.4–0.8 ${M}_{\odot }$). Kojima-1Lb is placed at the region where only a little has yet been surveyed by any methods due to the limitation of their sensitivity. Several planets have been discovered in the same region with the radial velocity technique (e.g., Mordasini et al. 2011; Astudillo-Defru et al. 2017), which, however, provides only a lower limit on their masses. On the other hand, the absolute mass of Kojima-1Lb is measured with an uncertainty of only 10%.

Zoom In Zoom Out Reset image size

The orbit of Kojima-1Lb was likely comparable to the snow line at its younger age, when the planet probably formed from a protoplanetary disk. We estimate that the snow-line location in the protoplanetary disk of Kojima-1L is ∼1.6 au by using the conventional formula of ${a}_{\mathrm{snow}}=2.7\times {M}_{* }/{M}_{\odot }$ au (e.g., Bennett et al. 2008; Sumi et al. 2010; Muraki et al. 2011), where M_* is the stellar mass. This mass–linear relation can be derived by assuming that the stellar luminosity is proportional to ${M}_{* }^{2}$ and the protoplanetary disk is optically thin (Bennett et al. 2008). Under this simple assumption, the present location of Kojima-1Lb is comparable to or slightly inner from the snow-line location of its youth, as shown in Figure 11(b).

More realistically, the snow-line distance is a function of age due to the evolution of the protoplanetary disk and stellar luminosity (e.g., Kennedy et al. 2006; Kennedy & Kenyon 2008). In Figure 12, we compare the orbit of Kojima-1Lb with a theoretical prediction of the time evolution of the snow-line location at the midplane of a young disk around a 0.6 ${M}_{\odot }$ star by Kennedy & Kenyon (2008; extracted from Figure 1 of their paper). The model assumes stellar irradiation and viscous accretion as the sources of disk heating. According to this model, the snow-line distance monotonically decreases with time, crossing the current planet location at an age of ${2.2}_{-1.6}^{+1.7}$ Myr. This timescale is comparable to or shorter than the typical disk lifetime of low-mass stars of a few tens of Myr (e.g., Luhman & Mamajek 2012; Ribas et al. 2015), indicating that the current location of Kojima-1Lb could have experienced a period when it was outside the snow line while disk gas remained.

Zoom In Zoom Out Reset image size

According to the core accretion theories, it is difficult to form a planet as massive as Kojima-1Lb (20 ± 2 ${M}_{\oplus }$) inside the snow line because of the lack of materials (e.g., Ida & Lin 2005; Kennedy et al. 2006), unless the surface density of solid materials in the disk's inner region is substantially high (e.g., Hansen & Murray 2012; Ogihara et al. 2015). On the other hand, in situ formation of Kojima-1Lb would be possible during the period when the snow line was inside the orbit of Kojima-1Lb and the disk gas still remained. Solid materials are thought to be abundant around the snow line (e.g., Kokubo & Ida 2002; Dra̧żkowska & Alibert 2017), which would allow the protoplanet of Kojima-1Lb to reach a mass of several ${M}_{\oplus }$ and start to accrete the surrounding gas. Several population-synthesis studies including type I migration also predict efficient formation of Neptune-mass planets near the snow line (e.g., Ida & Lin 2005; Mordasini et al. 2009), while the recent result of microlensing surveys has required some modifications of these predictions, at least for the region outside a few times the snow line (Suzuki et al. 2018). Although it is not possible to identify the exact formation process of this specific planet, given the precise mass determination of Kojima-1Lb, this planet could be an important example toward understanding the planetary formation processes around the snow line.

It is interesting to consider the detection efficiency of the planetary signal in Kojima-1, as the sensitivity to the planet in this event could be different from typical microlensing events toward the Galactic bulge.

Assuming that the actual planet signal is absent, we calculate the detection efficiency by following the method of Rhie et al. (2000). In this calculation, we use not only the data sets that are used for the light-curve fitting but also all of the other data sets listed in Table 1, except for the SL data set that was identified to have systematics. On the other hand, we eliminate all data points after 2018 January 1 (HJD–2,450,000 = 8120), because we would have terminated our photometric follow-up campaign by the end of 2017 if the planetary signal was not detected. As the threshold of signal detection, we adopt ${\rm{\Delta }}{\chi }^{2}=100$ following Suzuki et al. (2016), where ${\rm{\Delta }}{\chi }^{2}$ is the ${\chi }^{2}$ difference between planetary and nonplanetary (single-lens) models. At first, the detection efficiency is computed as a function of $(\mathrm{log}\,s,\mathrm{log}\,q)$. Next, we transform it to the physical parameter space, $(\mathrm{log}\,{a}_{\mathrm{proj}},\mathrm{log}\,{M}_{L2})$ (Dominik 2006), where we use the well-constrained probability distribution function of ${\theta }_{E}$ and M_L1 instead of the Bayesian approach using a Galactic model. The detection efficiency $\epsilon (\mathrm{log}\,{a}_{\mathrm{proj}},\mathrm{log}\,{M}_{L2})$ is further converted to $\epsilon (\mathrm{log}\,{a}_{3{\rm{D}}},\mathrm{log}\,{M}_{L2})$ with the assumption that the planet has a circular orbit and random orientation.

The calculated detection efficiency is plotted by contours in Figure 11(a). We also calculate the detection efficiency as a function of $\mathrm{log}({a}_{3{\rm{D}}}/{a}_{\mathrm{snow}})$ and $\mathrm{log}\,{M}_{L2}$, where ${a}_{\mathrm{snow}}=2.7\,\times ({M}_{* }/{M}_{\odot })\mathrm{au}$, as shown in Figure 11(b). The planet sensitivity of Kojima-1 has its peak around 1–1.4 au, or 0.7–1.0 times the snow-line distance. This region is a few times interior to the region where the majority of microlensing planets have been discovered, reflected by the fact that the distance to the source star of Kojima-1 is ∼10 times closer to us than those of the other microlensing events.

On the other hand, the detection efficiency of Kojima-1Lb is calculated to be only ∼35%. Here we remind the reader that the Kojima-1 event was not discovered by a systematic microlensing survey but was unexpectedly discovered during a nova search conducted by an amateur astronomer. Only one such event was previously known (the so-called Tago event; Fukui et al. 2007; Gaudi et al. 2008), but in that case, no planetary signal was detected. Therefore, although it is too early to argue statistically, the discovery of this low detection efficiency planet may imply that Neptunes are common rather than rare in this orbital region. This result is consistent with the recent findings with transit and radial velocity techniques that Neptunes are at least as common as (Kawahara & Masuda 2019) or more common than (Herman et al. 2019; Tuomi et al. 2019) Jupiters at large orbits comparable to the snow line.

Unlike many of the other microlensing planetary systems, Kojima-1L offers valuable opportunities to follow up in various ways thanks to its closeness to the Earth. First, the geocentric source-lens relative proper motion is estimated to be ${\mu }_{\mathrm{geo}}=25.34\pm 0.44$ mas yr⁻¹, enabling us to spatially separate the source and lens stars in ∼2 yr from the event using ground-based AO instruments (e.g., Keck/NIRC2) or space-based telescopes (e.g., Hubble Space Telescope). By resolving the two stars, one can confirm the relative proper motion (including its direction) and the brightness of the host star in an independent way (e.g., Batista et al. 2015; Bennett et al. 2015; Bhattacharya et al. 2018).

Second, the host star is as bright as K_s = 13.7, which is the brightest among all known microlensing planetary systems followed by OGLE-2018-BLG-0740L (Han et al. 2019), allowing spectroscopic characterizations of the host star. Low- or mid-resolution spectroscopy in the near-infrared is feasible with a >4 m class telescope, ideally with an AO instrument to reduce the contamination flux from the background source star. Such an observation will provide fundamental spectroscopic information on the host star, such as temperature, metallicity, and kinematics in the Galaxy. Furthermore, it is possible to search for additional inner and/or more massive planets with the radial velocity technique using an 8 m class telescope equipped with an AO-guided, near-infrared, high-dispersion spectrograph, such as Subaru/IRD. Knowing planetary multiplicity is of particular importance in understanding the formation and dynamical evolution of this planetary system. Finally, Kojima-1Lb would induce a radial velocity on the host star with an amplitude of ∼2.2 sin i ms⁻¹ and a period of ∼1.5 yr assuming a circular orbit, where i is orbital inclination. This signal will be measurable in the era of extremely large telescopes (ELTs), offering a valuable opportunity to confirm the mass and refine the orbit of this snow-line Neptune.