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

A. Fukui, D. Suzuki, N. Koshimoto, E. Bachelet, T. Vanmunster, D. Storey, H. Maehara, K. Yanagisawa, T. Yamada, A. Yonehara, et al.

► To cite this version:

A. Fukui, D. Suzuki, N. Koshimoto, E. Bachelet, T. Vanmunster, et al. Kojima-1Lb Is a Mildly Cold Neptune around the Brightest Microlensing Host Star. The Astronomical Journal, 2019, 158 (5), pp.206. 10.3847/1538-3881/ab487f. hal-02404798

HAL Id: hal-02404798

https://hal.science/hal-02404798v1

Submitted on 4 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

A. Fukui<sup>1,2,●</sup>, D. Suzuki<sup>3</sup>, N. Koshimoto<sup>4,5,6</sup>, E. Bachelet<sup>7,●</sup>, T. Vanmunster<sup>8</sup>, D. Storey<sup>9</sup>, H. Maehara<sup>10</sup>, K. Yanagisawa<sup>11</sup>,

T. Yamada<sup>12</sup>, A. Yonehara<sup>12</sup>, T. Hirano<sup>13</sup>, D. P. Bennett<sup>6,14</sup>, V. Bozza<sup>15,16</sup>, D. Mawet<sup>17,18,●</sup>, M. T. Penny<sup>19,●</sup>, S. Awiphan<sup>20</sup>,

A. Oksanen<sup>9,21</sup>, T. M. Heintz<sup>22,23</sup>, T. E. Oberst<sup>23</sup>, V. J. S. Béjar<sup>2,24</sup>, N. Casasayas-Barris<sup>2,24</sup>, G. Chen<sup>2,25</sup>, N. Crouzet<sup>2,24,26</sup>,

D. Hidalgo<sup>2,24</sup>, P. Klagyikiv<sup>2,24</sup>, F. Murgas<sup>2,24</sup>, N. Narita<sup>2,5,27,28,●</sup>, E. Palle<sup>2,24</sup>, H. Parviainen<sup>2,24</sup>, N. Watanabe<sup>5,29</sup>,

N. Kusakabe<sup>5,27</sup>, M. Mori<sup>4</sup>, Y. Terada<sup>4</sup>, J. P. de Leon<sup>4</sup>, A. Hernandez<sup>2,24</sup>, R. Luque<sup>2,24</sup>, M. Monelli<sup>2,24</sup>, P. Montañes-Rodriguez<sup>2,24</sup>, J. Prieto-Arranz<sup>2,24</sup>, K. L. Murata<sup>30</sup>, S. Shugarov<sup>31,32</sup>, Y. Kubota<sup>12</sup>, C. Otsuki<sup>12</sup>, A. Shionoya<sup>12</sup>, T. Nishiumi<sup>5,12</sup>, A. Nishide<sup>12</sup>,

M. Fukagawa<sup>5,●</sup>, K. Onodera<sup>3,33,34,35</sup>, S. Villanueva, Jr.<sup>36</sup>, R. A. Street<sup>7</sup>, Y. Tsapras<sup>37,●</sup>, M. Hundertmark<sup>37,●</sup>, M. Kuzuhara<sup>27</sup>, M. Fujita<sup>13</sup>, C. Beichman<sup>18,38,39</sup>, J.-P. Beaulieu<sup>40,41</sup>, R. Alonso<sup>2,24</sup>, D. E. Reichart<sup>42,●</sup>, N. Kawai<sup>30</sup>, and M. Tamura<sup>4,5,27,●</sup>

<sup>1</sup>Department of Earth and Planetary Science, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

<sup>2</sup>Instituto de Astrofísica de Canarias, Vía Láctea s/n, E-38205 La Laguna, Tenerife, Spain

<sup>3</sup>Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency (JAXA), 3-1-1 Yoshinodai, Chuo, Sagamihara, Kanagawa 252-5210, Japan

<sup>4</sup>Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

<sup>5</sup>National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

<sup>6</sup>Code 667, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

<sup>7</sup>Las Cumbres Observatory, 6740 Cortona Drive, Suite 102, Goleta, CA 93117, USA

<sup>8</sup>Center for Backyard Astrophysics Belgium, Walhostraat 1A, B-3401 Landen, Belgium

<sup>9</sup>American Association of Variable Star Observers, 49 Bay State Road, Cambridge, MA 02138, USA

<sup>10</sup>Okayama Branch Office, Subaru Telescope, National Astronomical Observatory of Japan, NINS, Kamogata, Asakuchi, Okayama 719-0232, Japan

<sup>11</sup>Division of Optical and Infrared Astronomy, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka-shi, Tokyo 181-8588, Japan

<sup>12</sup>Department of Astrophysics and Atmospheric Sciences, Faculty of Science, Kyoto Sangyo University, 603-8555 Kyoto, Japan

<sup>13</sup>Department of Earth and Planetary Sciences, Tokyo Institute of Technology, 2-12-1 Okayama, Meguro-ku, Tokyo 152-8551, Japan

<sup>14</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA

<sup>15</sup>Dipartimento di Fisica "E.R. Caianiello," Università di Salerno, Via Giovanni Paolo II 132, I-84084, Fisciano, Italy

<sup>16</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Napoli, Via Cintia, I-80126 Napoli, Italy

<sup>17</sup>Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

<sup>18</sup>Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA

<sup>19</sup>Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, OH 43210, USA

<sup>20</sup>National Astronomical Research Institute of Thailand, 260, Moo 4, T. Donkaew, A. Mae Rim, Chiang Mai, 50180, Thailand

<sup>21</sup>Hankasalmi Observatory, Hankasalmi, Finland

<sup>22</sup>Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA

<sup>23</sup>Department of Physics, Westminster College, New Wilmington, PA 16172, USA

<sup>24</sup>Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38206 La Laguna, Tenerife, Spain

<sup>25</sup>Key Laboratory of Laboratory of Astrophysics Sciences, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008, People's Republic of China

<sup>26</sup>European Space Agency, European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

<sup>27</sup>Astrobiology Center, National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

<sup>28</sup>JST, PRESTO, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

<sup>29</sup>SOKENDAI (The Graduate University of Advanced Studies), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

<sup>30</sup>Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan

<sup>31</sup>Sternberg Astronomical Institute, Moscow State University, Moscow 119991, Russia

<sup>32</sup>Astronomical Institute of the Slovak Academy of Sciences, Tatranska Lomnica, 05960, Slovakia

<sup>33</sup>Department of Space and Astronautical Science, SOKENDAI (The Graduate University for Advanced Studies), 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210, Japan

<sup>34</sup>Institut de Physique du Globe de Paris, 1 Rue Jussieu, F-75005 Paris, France

<sup>35</sup>Université Paris Diderot, 5 Rue Thomas Mann, F-75013 Paris, France

<sup>36</sup>Kavli Institute for Astrophysics and Space Research, M.I.T., Cambridge, MA 02139, USA

<sup>37</sup>Zentrum für Astronomie der Universität Heidelberg, Astronomisches Rechen-Institut, Mönchhofstr. 12-14, D-69120 Heidelberg, Germany

<sup>38</sup>Division of Physics, Mathematics, and Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

<sup>39</sup>NASA Exoplanet Science Institute, 770 South Wilson Avenue, Pasadena, CA 911225, USA

<sup>40</sup>School of Physical Sciences, University of Tasmania, Private Bag 37 Hobart, Tasmania 7001, Australia

<sup>41</sup>Sorbonne Universites, CNRS, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France

<sup>42</sup>Department of Physics and Astronomy, University of North Carolina, CB #3255, Chapel Hill, NC 27599, USA

Received 2019 May 24; revised 2019 September 16; accepted 2019 September 25; published 2019 October 31

## Abstract

We report the analysis of additional multiband photometry and spectroscopy and new adaptive optics (AO) imaging of the nearby planetary microlensing event TCP J05074264+2447555 (Kojima-1), which was discovered toward the Galactic anticenter in 2017 (Nucita et al.). We confirm the planetary nature of the light-curve anomaly around the peak while finding no additional planetary feature in this event. We also confirm the presence of apparent blending flux and the absence of significant parallax signal reported in the literature. The AO image reveals no contaminating sources, making it most likely that the blending flux comes from the lens star. The measured multiband lens flux, combined with a constraint from the microlensing model, allows us to narrow down the previously unresolved mass and distance of the lens system. We find that the primary lens is a dwarf on the K/M boundary (0.581 ± 0.033 $M_\odot$) located at 505 ± 47 pc, and the companion (Kojima-1Lb) is a Neptune-mass planet (20.0 ± 2.0 $M_\oplus$) with a semimajor axis of 1.08$^{+0.62}_{-0.18}$ au. This orbit is a few times smaller than those of typical microlensing planets and is comparable to the snow-line location at young ages. We calculate that the a priori
detection probability of Kojima-1Lb is only ~35%, which may imply that Neptunes are common around the snow line, as recently suggested by the transit and radial velocity techniques. The host star is the brightest among the microlensing planetary systems ($K_s = 13.7$), offering a great opportunity to spectroscopically characterize this system, even with current facilities.

Unified Astronomy Thesaurus concepts: Gravitational microlensing (672); Exoplanet systems (484)

# 1. Introduction

According to core accretion theory, once a protoplanetary core reaches a critical mass of ~10 $M_⊕$ by accumulating planetesimals, the protoplanet starts to accrete the surrounding gas in a runaway fashion and quickly becomes a gas giant planet (e.g., Pollack et al. 1996). This process can most efficiently happen just outside the snow line, where the surface density of solid materials is enhanced by condensation of ices (e.g., Ida & Lin 2004). Because this process is basically controlled by the mass of the protoplanet, unveiling the planetary mass distribution around the snow line is crucial to understand the planetary formation processes. Recent microlensing surveys have revealed that Neptune-mass-ratio planets are the most abundant in the region several times outside the snow line (Suzuki et al. 2016; Udalski et al. 2018); however, little is known about the population of low-mass planets just around the snow line.

The microlensing technique is most sensitive to planets with an orbital separation close to the Einstein radius, which is defined by the radius of the ringed image produced when the lens and source stars are perfectly aligned. This size is expressed by

$$R_E = \sqrt{\frac{4G}{c^2} M_L D_S x(1-x)} \quad (1)$$

$$\approx 2.9 \text{ au} \left(\frac{M_L}{0.5M_\odot}\right)^{1/2} \left(\frac{D_S}{8 \text{ kpc}}\right)^{1/2} \left[\frac{x(1-x)}{0.25}\right]^{1/2}, \quad (2)$$

where $M_L$ is the mass of the lens star, $x = D_L/D_S$, and $D_L$ and $D_S$ are the distances to the lens and source stars, respectively. Assuming that the snow-line distance in a protoplanetary disk can be approximated by $a_{\text{snow}} \sim 2.7 \text{ au} \times M_*/M_⊕$, where $M_*$ is the stellar mass (Bennett et al. 2008), one can write the ratio of the Einstein radius to the median sky-projected distance of the randomly oriented snow-line orbit, $a_{\text{snow,⊥}} = 0.866a_{\text{snow}}$, as

$$\frac{R_E}{a_{\text{snow,⊥}}} \approx 2.4 \left(\frac{M_L}{0.5M_\odot}\right)^{-1/2} \left(\frac{D_S}{8\text{kpc}}\right)^{1/2} \left[\frac{x(1-x)}{0.25}\right]^{1/2}. \quad (3)$$

Thus, the Einstein radius of typical microlensing events toward the Galactic bulge ($M_L \sim 0.5 M_⊕$, $x \sim 0.5$, and $D_S \sim 8$ kpc), where dedicated microlensing surveys have been conducted, is a few times larger than the snow-line distance (see, e.g., Tsapras 2018, for a recent review of microlensing).

Because the Einstein radius is scaled by $\sqrt{D_S}$, the planet sensitivity region of microlensing coincides with the location of the snow line when the distance of the source is an order of magnitude closer than the distance to the Galactic bulge, i.e., $D_S \sim 1$ kpc. Although the event rate of such nearby-source microlensing events is expected to be small (~23 events yr$^{-1}$; Han 2008), they can provide a rare opportunity to find and

characterize planets just around the snow line. In addition, once such a nearby planetary microlensing event is discovered, it can be an invaluable system that allows spectroscopic follow-up, which is usually difficult for the events observed toward the Galactic bulge.

This is the case for the nearby microlensing event TCP J05074264+2447555$^{43}$ (hereafter Kojima-1$^{44}$), which was serendipitously discovered during a nova search conducted by an amateur astronomer, Mr. T. Kojima. On 2017 October 31 UT, he reported an unknown transient event on an $R = 13.6$ mag star toward the Taurus constellation,$^{45}$ and later, the microlensing nature of this event was confirmed by photometric and spectroscopic follow-up observations (Jayasinghe et al. 2017; Konyves-Toth et al. 2017; Maehara 2017; Sokolovsky 2017). Moreover, a planetary feature was detected near the peak of the event by the earliest photometric follow-up observations (Nucita et al. 2017).

Nucita et al. (2018) estimated that the distance to the source star is ~700–800 pc. They also fit their own and publicly available light curves with a binary-lens microlens model, finding that the mass ratio of the primary lens to its companion is $(1.1 \pm 0.1) \times 10^{-4}$; i.e., the companion is a planet. However, because of the degeneracy between the absolute mass and distance of the lens system, they estimated them using a stochastic technique based on a Galactic model such that the planetary mass is $9.2 \pm 6.6 M_⊕$, the host star’s mass is ~0.25 $M_⊕$, and the distance to the system is ~380 pc. On the other hand, Dong et al. (2019) measured the angular Einstein radius $\theta_E$ of this event by observing the separation of the two microlensed source star images using the VLTI/GRAVITY instrument. They confirmed that the $\theta_E$ value estimated by Nucita et al. (2018) is largely consistent with the value measured by VLTI, although they did not attempt to improve the physical parameters of the lens system using the improved $\theta_E$.

Reacting to the discovery of this remarkable event, we started follow-up observations by means of photometric monitoring, high- and low-resolution spectroscopy, and high-resolution imaging to obtain a better understanding of the lens system.

This paper is organized as follows. We describe our follow-up observations and reductions in Section 2 and light-curve modeling in Section 3. The properties of the source star and lens system are derived in Sections 4 and 5, respectively. We then discuss the possible formation scenario of the planet, detection efficiency of the planet, and capabilities of future

<footnote>$^{43}$ The equatorial and galactic coordinates of this object are $(\alpha, \delta)_{2000} = (05^h07^m42^s725, +24^s47^m56.33'37)$ and $(l, b)_{2000} = (178^\circ76, -9^\circ32)$, respectively.</footnote>

<footnote>$^{44}$ Note that Nucita et al. (2018) nicknamed this event Feynman-01 in honor of the observatory where the planetary feature was observed. In this paper, we call this event Kojima-1 in honor of Mr. Kojima as the first discoverer of this event. Conventionally, a planetary microlensing event is named after the group(s) that discovers the event itself, rather than the group(s) that detects the planetary feature.</footnote>

<footnote>$^{45}$ http://www.cbat.eps.harvard.edu/unconf/followups/J05074264+2447555.html</footnote>
<table><caption>**Table 1**
List of Photometric Data Sets</caption><thead><tr><th>Abbreviation<br>(Instrument)<sup>a,b</sup></th><th>Observatory</th><th>Telescope<br>Diameter<br>(m)</th><th>Field of View<br>(arcmin<sup>2</sup>)</th><th>Filter</th><th>Number of<br>Nights<sup>c</sup></th><th>Number of<br>Data<sup>c</sup></th><th>Median<br>Cadence<sup>c</sup><br>(minutes)</th><th>Median<br>Flux Error<sup>c</sup><br>(%)</th></tr></thead><tbody><tr><td colspan="9"><em>Data Sets Obtained or Rereduced in This Work</em></td></tr><tr><td>MuSCAT</td><td>NAOJ/Okayama</td><td>1.88</td><td>6.1 × 6.1</td><td>g</td><td>11</td><td>161</td><td>10.0</td><td>0.24</td></tr><tr><td></td><td></td><td></td><td></td><td>r</td><td>12</td><td>163</td><td>10.0</td><td>0.16</td></tr><tr><td></td><td></td><td></td><td></td><td>z<sub>s</sub></td><td>12</td><td>196</td><td>10.0</td><td>0.30</td></tr><tr><td>MuSCAT2</td><td>Teide Observatory</td><td>1.52</td><td>7.4 × 7.4</td><td>g</td><td>29</td><td>331</td><td>10.0</td><td>0.38</td></tr><tr><td></td><td></td><td></td><td></td><td>r</td><td>27</td><td>317</td><td>10.0</td><td>0.21</td></tr><tr><td></td><td></td><td></td><td></td><td>i</td><td>29</td><td>316</td><td>10.0</td><td>0.21</td></tr><tr><td></td><td></td><td></td><td></td><td>z<sub>s</sub></td><td>30</td><td>343</td><td>10.0</td><td>0.24</td></tr><tr><td>Araki</td><td>Koyama Astronomical Observatory</td><td>1.3</td><td>12.2 × 12.2</td><td>g</td><td>12</td><td>68</td><td>12.1</td><td>0.29</td></tr><tr><td></td><td></td><td></td><td></td><td>R<sub>c</sub></td><td>12</td><td>70</td><td>10.8</td><td>0.56</td></tr><tr><td>ISAS</td><td>JAXA/ISAS</td><td>1.3</td><td>5.4 × 5.4</td><td>I<sub>c</sub></td><td>8</td><td>175</td><td>10.1</td><td>0.67</td></tr><tr><td>OAOWFC</td><td>NAOJ/Okayama</td><td>0.91</td><td>28.6 × 28.6</td><td>K<sub>s</sub></td><td>43</td><td>202</td><td>56.0</td><td>1.95</td></tr><tr><td>CBABO</td><td>CBA Belgium Observatory</td><td>0.40</td><td>12.5 × 8.4</td><td>Clear</td><td>5</td><td>30</td><td>4.9</td><td>0.77</td></tr><tr><td>COAST</td><td>Teide Observatory</td><td>0.35</td><td>33 × 33</td><td>V</td><td>6</td><td>7</td><td>...</td><td>1.18</td></tr><tr><td>PROMPT-8</td><td>Cerro Tololo Inter-American Observatory</td><td>0.61</td><td>22.6 × 22.6</td><td>V</td><td>8</td><td>64</td><td>9.7</td><td>0.83</td></tr><tr><td></td><td></td><td></td><td></td><td>R<sub>c</sub></td><td>9</td><td>79</td><td>9.7</td><td>0.69</td></tr><tr><td></td><td></td><td></td><td></td><td>I<sub>c</sub></td><td>7</td><td>70</td><td>9.7</td><td>1.21</td></tr><tr><td>SL</td><td>AISAS in Stará Lesná</td><td>0.60</td><td>14.4 × 14.4</td><td>B</td><td>3</td><td>114</td><td>4.6</td><td>2.08</td></tr><tr><td></td><td></td><td></td><td></td><td>V</td><td>3</td><td>198</td><td>5.2</td><td>1.18</td></tr><tr><td></td><td></td><td></td><td></td><td>R<sub>c</sub></td><td>3</td><td>121</td><td>4.8</td><td>1.63</td></tr><tr><td></td><td></td><td></td><td></td><td>I<sub>c</sub></td><td>3</td><td>177</td><td>5.2</td><td>1.37</td></tr><tr><td>MITSuME</td><td>NAOJ/Okayama</td><td>0.50</td><td>26 × 26</td><td>I<sub>c</sub></td><td>28</td><td>239</td><td>13.3</td><td>1.06</td></tr><tr><td>DEMONEXT</td><td>Winer Observatory</td><td>0.50</td><td>30.7 × 30.7</td><td>I<sub>c</sub></td><td>20</td><td>420</td><td>10.5</td><td>2.73</td></tr><tr><td>OAR</td><td>Hankasalmi Observatory</td><td>0.40</td><td>25 × 25</td><td>V</td><td>4</td><td>39</td><td>6.4</td><td>0.68</td></tr><tr><td>WCO</td><td>Westminster College Observatory</td><td>0.35</td><td>24 × 16</td><td>CBB</td><td>5</td><td>129</td><td>9.8</td><td>0.18</td></tr><tr><td colspan="9"><em>Public or Published Data Sets</em></td></tr><tr><td>FO</td><td>R.P. Feynman Observatory</td><td>0.30</td><td>27.0 × 21.6</td><td>V</td><td>5</td><td>54</td><td>8.8</td><td>0.59</td></tr><tr><td>ASAS-SN</td><td>Haleakala Observatory</td><td>0.14</td><td>273 × 273</td><td>V</td><td>44</td><td>146</td><td>...</td><td>2.27</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The data sets used in the light-curve fitting are shown in bold.

<sup>b</sup> References to the instruments are as follows. MuSCAT: Narita et al. (2015); MuSCAT2: Narita et al. (2019); OAOWFC: Yanagisawa et al. (2016); MISTuME: Kotani et al. (2005), Yanagisawa et al. (2010); DEMONEXT: Villanueva et al. (2018).

<sup>c</sup> The values for the data after removing outliers and binning time series are reported.

follow-up observations of the planetary system in Section 6. We summarize the paper in Section 7.

## 2. Observations

### 2.1. Photometric Monitoring

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,
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (First panel) Light curves of Kojima-1. Colored (including black) and light gray points are the data used for the light-curve fitting and only for the calculation of detection efficiency, respectively. Color legends are shown on the left-hand side. The best-fit microlensing model is indicated by a blue solid line. The times when the two LCO spectra were taken are indicated by arrows. (Second panel) Residuals from the best-fit model. (Third panel) Zoomed light curves around the peak. The time when the HIDES spectrum was obtained is indicated by an arrow. (Fourth panel) Residuals for the zoomed light curves.</figcaption></figure>

ISAS, OAOWFC, CBABO, COAST, SL, and MITSuME; IRAF/APPHOT<fnref n="46" /> 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<fnref n="47" /> (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.

## 2.2. High-resolution Spectroscopy

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.

## 2.3. Low-resolution Spectroscopy

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.<fnref n="48" /> The spectral range is about 3200–10000 Å. Each spectrum was taken with 1000 s exposure with the 1′′2 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.<fnref n="49" />

## 2.4. High-resolution Imaging

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*<sub>s</sub> 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 0′′06. 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*<sub>s</sub> = 21 within the image.

<footnote n="46">⁴⁶ IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.</footnote>

<footnote n="47">⁴⁷ http://www2.iap.fr/users/alard/package.html</footnote>

<footnote n="48">⁴⁸ More details on the LCO instruments and telescope are available here: https://lco.global/observatory/.</footnote>

<footnote n="49">⁴⁹ https://github.com/svalenti/FLOYDS_pipe</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (Left) The K<sub>x</sub>-band AO image of the Kojima-1 object obtained with Keck/NIRC2. (Right) A 5σ contrast curve as a function of the distance from the centroid of the object.</figcaption></figure>

### 3. Light-curve Modeling

#### 3.1. Model Description

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<sub>0</sub>; the Einstein radius crossing time, t<sub>E</sub>; the source-lens angular separation at time t<sub>0</sub> in units of the angular Einstein radius (θ<sub>E</sub>), u<sub>0</sub>; the mass ratio of the binary components, q; the sky-projected separation of the binary components in units of θ<sub>E</sub>, s; the angle between the source trajectory and the binary-lens axis, α; the angular source radius in units of θ<sub>E</sub>, ρ; and the microlens parallax vector, π<sub>E</sub>. Here the direction of π<sub>E</sub> is the same as the direction of the source's proper motion relative to the lens, and the length of π<sub>E</sub>, π<sub>E</sub> = √(π<sub>E,N</sub><sup>2</sup> + π<sub>E,E</sub><sup>2</sup>), is equal to the ratio of 1 au to the projected Einstein radius onto the observer plane, where π<sub>E,N</sub> and π<sub>E,E</sub> are the north and east components of π<sub>E</sub>, respectively. The limb-darkening effect of the source star is modeled by the formula I(θ) = I(0)[1 - u<sub>x</sub>(1 - cos θ)], 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<sub>x</sub> 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<sub>i</sub>(t) = A(t) × F<sub>s,i</sub> + F<sub>b,i</sub>, where F<sub>s,i</sub> and F<sub>b,i</sub> 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.

#### 3.2. Error Normalization

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

$$ \sigma_i' = k \sqrt{\sigma_i^2 + e_{\min}^2}, \qquad (4) $$

where σ<sub>i</sub> is the initial uncertainty of the ith data point in magnitude, and k and e<sub>min</sub> are coefficients for each data set. Here the term e<sub>min</sub> represents systematic errors that dominate when the flux is significantly increased. The k and e<sub>min</sub> values are adjusted so that the cumulative χ<sup>2</sup> distribution for the best-fit binary-lens model including the parallax effect sorted by magnitude is close to linear and χ<sub>red</sub><sup>2</sup> 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 π<sub>E</sub>, 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.

#### 3.3. Data Sets and Fitting Codes

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 Super novae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017; data are extracted from their website<sup>50</sup> 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<sub>s</sub> and F<sub>b</sub> 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,<sup>51</sup> 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

<footnote n="50"><sup>50</sup> https://asas-sn.osu.edu</footnote>

<footnote n="51"><sup>51</sup> Although we found no clear evidence for the cause of this systematics, the stellar positions on the detector moved by >50 pixels during the observations, which might cause systematics on the photometry at some level.</footnote>
<table><caption>**Table 2**

Best-fit Parameter Values of Binary-lens Microlensing Models</caption><thead><tr><th>Parameter<sup>a</sup></th><th>Unit</th><th>Nucita et al. (2018)</th><th>Static</th><th>Parallax</th><th>Parallax<br>w/θ<sub>E</sub> Prior</th><th>Parallax<br>w/θ<sub>E</sub>, Φ<sub>π</sub> Priors</th></tr></thead><tbody><tr><td>t<sub>0</sub></td><td>HJD<br>-2,458,058</td><td>0.75 ± 0.01</td><td>0.7353 ± 0.0076</td><td>0.7395 ± 0.0073</td><td>0.7396 ± 0.0073</td><td>0.7403 ± 0.0074</td></tr><tr><td>t<sub>E</sub></td><td>days</td><td>26.4 ± 0.9</td><td>27.44 ± 0.07</td><td>27.19 ± 0.15</td><td>27.18 ± 0.14</td><td>27.25 ± 0.09</td></tr><tr><td>u<sub>0</sub></td><td>10<sup>-2</sup></td><td>9.3 ± 0.1<sup>b</sup></td><td>8.858<sup>+0.031</sup><sub>-0.034</sub></td><td>8.925 ± 0.043</td><td>8.927 ± 0.042</td><td>8.935 ± 0.038</td></tr><tr><td>q</td><td>10<sup>-4</sup></td><td>1.1 ± 0.1</td><td>1.058<sup>+0.068</sup><sub>-0.074</sub></td><td>1.075<sup>+0.066</sup><sub>-0.073</sub></td><td>1.031<sup>+0.078</sup><sub>-0.084</sub></td><td>1.027<sup>+0.078</sup><sub>-0.084</sub></td></tr><tr><td>s (model a)</td><td></td><td>0.935 ± 0.004</td><td>0.9207<sup>+0.0045</sup><sub>-0.0040</sub></td><td>0.9204<sup>+0.0040</sup><sub>-0.0038</sub></td><td>0.9263 ± 0.0018</td><td>0.9264 ± 0.0018</td></tr><tr><td>s (model b)</td><td></td><td>0.975 ± 0.004</td><td>0.9944<sup>+0.0041</sup><sub>-0.0046</sub></td><td>0.9941<sup>+0.0042</sup><sub>-0.0045</sub></td><td>0.9874 ± 0.0018</td><td>0.9873 ± 0.0018</td></tr><tr><td>α</td><td>rad</td><td>4.767 ± 0.007<sup>b</sup></td><td>4.7594 ± 0.0030</td><td>4.7610 ± 0.0030</td><td>4.7604 ± 0.0028</td><td>4.7604 ± 0.0028</td></tr><tr><td>ρ</td><td>10<sup>-3</sup></td><td>6.0 ± 0.8</td><td>3.2<sup>+0.9</sup><sub>-1.3</sub></td><td>3.2<sup>+0.9</sup><sub>-1.3</sub></td><td>4.568 ± 0.070</td><td>4.567 ± 0.071</td></tr><tr><td>π<sub>E,E</sub></td><td></td><td>...</td><td>...</td><td>0.071<sup>+0.072</sup><sub>-0.064</sub></td><td>0.0693<sup>+0.070</sup><sub>-0.063</sub></td><td>0.143<sup>+0.061</sup><sub>-0.053</sub></td></tr><tr><td>π<sub>E,N</sub></td><td></td><td>...</td><td>...</td><td>0.17 ± 0.45</td><td>0.19 ± 0.45</td><td>-0.33<sup>+0.12</sup><sub>-0.14</sub></td></tr><tr><td>χ<sup>2</sup><sub>min</sub> / dof</td><td></td><td>...</td><td>2557.4/2578</td><td>2543.0/2576</td><td>2546.5/2577</td><td>2550.7/2578</td></tr><tr><td>π<sub>E</sub></td><td></td><td>...</td><td>...</td><td>0.34<sup>+0.34c</sup><sub>-0.20</sub></td><td>0.35<sup>+0.34c</sup><sub>-0.20</sub></td><td>0.36<sup>+0.16</sup><sub>-0.13</sub></td></tr><tr><td>θ<sub>*</sub></td><td>μas</td><td>...</td><td>8.59 ± 0.06</td><td>8.65 ± 0.06</td><td>8.63 ± 0.06</td><td>8.63 ± 0.06</td></tr><tr><td>θ<sub>E</sub></td><td>mas</td><td>1.45 ± 0.25</td><td>2.63<sup>+1.77</sup><sub>-0.58</sub></td><td>2.68<sup>+1.87</sup><sub>-0.59</sub></td><td>1.890 ± 0.032</td><td>1.890 ± 0.032</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The 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.

<sup>b</sup> For ease of comparison, we multiply the *u*<sub>0</sub> and increment *α* reported in the literature by −1 and *π*, respectively. The geometry is identical to this transformation.

<sup>c</sup> Because $\pi_{E,E}$ and $\pi_{E,N}$ take both positive and negative values, the median value of $\pi_E$ does not coincide with $\sqrt{\langle\pi_{E,E}\rangle^2 + \langle\pi_{E,N}\rangle^2}$, where $\langle\pi_{E,E}\rangle$ and $\langle\pi_{E,N}\rangle$ are the median values of $\pi_{E,E}$ and $\pi_{E,N}$, respectively.

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<sup>52</sup> (Bozza et al. 2018) for sanity check.

### 3.4. Static Model

We first fit the light curves with a binary-lens model without the microlens parallax effect (static model), fixing $\pi_{EE}$ and $\pi_{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\sigma$ 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 $\chi^2$ 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.

Our derived values are consistent with those of Nucita et al. (2018) within $2\sigma$ for all parameters except for $u_0$, *s*, and $\rho$, 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. Parallax Model

#### 3.5.1. Without Informative Prior

To search for a signal of the parallax effect, we fit the light curves letting $\pi_{EE}$ and $\pi_{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 $\pi_E$ value of $0.34^{+0.34}_{-0.20}$. However, the $\chi^2$ 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 \ln N_{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 $\chi^2$ between the two models as a function of time. This plot indicates that most of the $\chi^2$ 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 $\pi_E$ (Section 3.5.3) and constrain the direction of $\pi_E$ (Section 3.5.4).

The result that a significant parallax signal is absent is consistent with the result of Dong et al. (2019), who also did

<footnote n="52">52 http://www.fisica.unisa.it/GravitationAstrophysics/RTModel.htm</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (Top) Difference of best-fit model magnitudes between the parallax and static models for individual data sets, where the color codes are the same as in Figure 1. (Bottom) Difference of cumulative $\chi^2$ between the parallax and static models for individual data sets (thin colored lines) and all data sets (thick gray line), where negative means that the parallax model is preferred. The color codes are the same as in Figure 1.</figcaption></figure>

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<fnref n="53" /> 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 $\theta_E$

From the light-curve fitting with the parallax model, $\rho$ is measured to be $3.2_{-1.3}^{+0.9} \times 10^{-3}$. This $\rho$ value allows the derivation of the angular Einstein radius $\theta_E$ via the relation of $\theta_E \equiv \theta_*/\rho$, where $\theta_*$ is the angular radius of the source star. The $\theta_*$ value is estimated to be $8.65 \pm 0.06$ µas using the procedure described in Section 4.4, which leads to $\theta_E = 2.7_{-0.6}^{+1.9}$ mas. On the other hand, the $\theta_E$ of the same event was independently and much more precisely determined to be $1.883 \pm 0.014$ mas (in the case of a luminous lens) by Dong et al. (2019) by spatially resolving the two microlensed images during the event. This information can be used to further constrain $\rho$ and some other parameters that are correlated with $\rho$ (in particular, $s$).

Using $\theta_E = 1.883 \pm 0.014$ mas (in the form of $\rho = \theta_*/\theta_E$) as an informative prior, we iteratively fit the light curves refining $\theta_*$ through the process described in Section 4.4. The improved parameter values are appended to Table 2, in which notable improvements can be seen in $\rho$, $s$, and $\theta_E$. On the other hand, the $\theta_E$ prior has not changed the significance of the parallax signal.

### 3.5.3. Upper Limit on $\pi_E$

From the VLTI observation, Dong et al. (2019) also constrained the direction of $\pi_E$ ($\Phi_\pi$) into two directions, $193^\circ.5 \pm 0^\circ.4$ and $156^\circ.7 \pm 0^\circ.4$ from north to east (for the luminous-lens model). To put an upper limit on $\pi_E$ utilizing the prior information of $\Phi_\pi$, we draw $\chi^2$ maps on a grid of $\pi_{E,E}$ and

$\pi_{E,N}$. We grid $\pi_{E,E}$ and $\pi_{E,N}$ by a grid size of 0.1 in the ranges of $-0.7 \le \pi_{E,E} < 0.7$ and $-1.5 \le \pi_{E,N} < 1.5$ and fit the light curves using the $\theta_E$ prior while fixing $\pi_{E,E}$ and $\pi_{E,N}$ at each grid-point value. In the left panel of Figure 4, we show $\Delta\chi^2$ maps on the $\pi_{E,E}-\pi_{E,N}$ plane calculated from all data sets, where $\Delta\chi^2$ is the difference of $\chi^2$ between each grid point and $(\pi_{E,E}, \pi_{E,N}) = (0, 0)$. The minimum-$\chi^2$ (darkest red) region is not coincident with the two solutions of $\Phi_\pi$ (indicated by cyan lines), probably due to the systematics in the light curves discussed before. Note that the reason why the negative $\Delta\chi^2$ region is elongated almost along the $\pi_{E,N}$ direction (only $\pi_{E,E}$ is well constrained) is that the direction of Earth’s acceleration is almost parallel to the direction of $\pi_{E,E}$.<fnref n="54" /> On the other hand, the right panel of the same figure shows a $\Delta\chi^2$ map that is calculated only using the $\chi^2$ 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-$\chi^2$ region is not localized, the intersection between the $\Phi_\pi$ solutions and some $\Delta\chi^2$ contour can still be used to put an upper limit on $\pi_E$. The contour of $\Delta\chi^2 = 9$ (white) intersects with the $\Phi_\pi \sim 156^\circ.7$ and $193^\circ.5$ lines (cyan) at the grid points that correspond to $\pi_E = 1.1$ and 0.5, respectively. We conservatively adopt 1.1 as a $3\sigma$ upper limit on $\pi_E$.

### 3.5.4. On the Direction of $\pi_E$

As will be discussed in Section 5.2.2, under the condition of $\pi_E < 1.1$, it is most likely that the blending flux detected in the light curves comes from the lens star independently on the $\Phi_\pi$ 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 $\theta_E$, predicts the $\pi_E$ value using the relation

$$ \pi_E = \frac{\theta_E}{\kappa M_L}, \qquad (5) $$

where $\kappa \equiv 4G/c^2$, $G$ is the gravitational constant, and $c$ is the speed of light. This gives $\pi_E = 0.39_{-0.03}^{+0.04}$, which is indicated by magenta solid (median) and dotted (1$\sigma$ boundary) contours in Figure 4. In the $\Delta\chi^2$ map for all data sets (left panel of Figure 4), the $\Delta\chi^2$ value at the grid point that satisfies both $\pi_E \sim 0.39$ and $\Phi_\pi \sim 156^\circ.7$ is $-16$, which is smaller than the counterpart that satisfies both $\pi_E \sim 0.39$ and $\Phi_\pi \sim 193^\circ.5$ by 40. This $\chi^2$ difference nominally rules out the $\Phi_\pi = 193^\circ.5$ solution.

This outcome, however, could be affected by systematics in the light curves. To test this possibility, we also check the $\Delta\chi^2$ map calculated using only the $\chi^2$ values from the ASAS-SN data set (right panel of Figure 4). We find that the $\Phi_\pi = 156^\circ.7$ solution is preferred over the other solution with a $\chi^2$ improvement of ~5, which, although marginal, supports the outcome obtained from all data sets.

Considering the above evidence, we adopt the $\Phi_\pi = 156^\circ.7$ solution for further analysis. To derive the final posteriors of the parameters, taking into account correlations between the parallax parameters ($\pi_{E,E}$ and $\pi_{E,N}$) and others and using all of the informative prior information, we rerun the MCMC analysis, letting $\pi_{E,E}$ and $\pi_{E,N}$ be free and imposing priors on $\theta_E$ and $\Phi_\pi$ with Gaussian distributions of $\theta_E = 1.883 \pm 0.014$ mas and

<footnote n="53">53 Dong et al. (2019) erroneously stated it to be 3 weeks.</footnote>

<footnote n="54">54 The ecliptic coordinate of the event is $(\beta, \lambda) = (78^\circ, 1^\circ 9')$, which is close to $(90^\circ, 0^\circ)$ where the direction of Earth’s acceleration is parallel to east-west.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (Left) The $\Delta\chi^2$ map for $\pi_{E,E}$ and $\pi_{E,E}$, where $\Delta\chi^2$ is the $\chi^2$ difference between each grid point and $(\pi_{E,E}, \pi_{E,N}) = (0, 0)$, calculated using all data sets. The two $\Phi_\pi$ solutions derived from the VLTI observation by Dong et al. (2019) are indicated by cyan lines. The magenta solid and dotted circles correspond to the contours of $\pi_E = 0.39_{-0.03}^{+0.04}$, which are expected from the lens flux (see text for details). (Right) Same as the left panel but calculated only using the $\chi^2$ of the ASAS-SN data set. The white solid lines are the contour for $\Delta\chi^2 = 9$. We estimate the $3\sigma$ upper limit of $\pi_E$ to be 1.1 from the intersection between the white and cyan lines.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Caustic (red) and source trajectory (gray) of the two degenerated microlensing models *a* (top) and *b* (bottom). The time ticks are given by small gray circles. The blue circle represents the source size and position at time $t = t_0$.</figcaption></figure>

$\Phi_\pi = 156.7 \pm 0.4$. 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 $\Phi_\pi$ 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 $\Phi_\pi$ value can be confirmed in the future by

directly measuring the lens-source relative position from high spatial resolution images.

## 4. Properties of the Source Star

In this section, we will derive the properties of the source star, in particular, the source's angular radius $\theta_*$ and the distance to the source star $D_S$, the former of which is tied to $\theta_E$ by the relation of $\theta_E = \theta_*/\rho$. We measure these values from the brightness of the source star derived from the light-curve fitting with the aid of the spectroscopic information and the extinction from the Gaia Data Release 2 (DR2).

### 4.1. High-resolution Spectrum

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_{\text{eff}} = 6303 \pm 110$ K, $R_S = 1.56 \pm 0.25$ R$_\odot$, and [Fe/H] = −0.11 ± 0.08, respectively. This result indicates that the source star is a main-sequence late-F dwarf.

### 4.2. Low-resolution Spectrum

The two LCO spectra were taken at magnifications of $A_1 = 8.34$ and $A_2 = 1.04$, with which the flux contamination from the lens star, in particular for the wavelength of $\gtrsim 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,\lambda}$ and $f_{2,\lambda}$ are the fluxes at the wavelength $\lambda$ in the first- and second-epoch
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Low-resolution spectrum of the source star extracted and extinction-corrected from the LCO spectra (green), along with empirical spectral templates of F4V, F5V, and F6V stars from Kesseli et al. (2017) (black, top to bottom).</figcaption></figure>

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.

### 4.3. Extinction Estimated from the Gaia DR2

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 ($\pi$) 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 $\pi$ and $A_G$, and have $\pi > 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\sigma$ error (standard deviation divided by the square root of the number of data points) for each bin are calculated, where the median $1\sigma$ error is ~0.10. The binned data are then fitted with a fourth-order polynomial function of the distance, which gives

$$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, \qquad (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$.

### 4.4. Distance and Angular Radius

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

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Extinction in the Gaia (left-hand axis; $A_G$) or visible (right-hand axis; $A_V$) band as a function of distance for stars within 30' in radius from the source position extracted from the Gaia DR2. Blue dots are the data for individual stars, and black squares are the binned values with a bin size of 50 pc, where the error bars represent the standard deviation divided by the square root of the number of data points. The red curve indicates the best-fit, fourth-order polynomial function.</figcaption></figure>

<table><caption>**Table 3**
Properties of the Source Star</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>g'</td><td>mag</td><td>14.559 ± 0.010</td></tr><tr><td>V</td><td>mag</td><td>14.151 ± 0.005</td></tr><tr><td>r'</td><td>mag</td><td>13.847 ± 0.008</td></tr><tr><td>i'</td><td>mag</td><td>13.556 ± 0.010</td></tr><tr><td>z'</td><td>mag</td><td>13.376 ± 0.009</td></tr><tr><td>K<sub>s</sub></td><td>mag</td><td>11.990 ± 0.012</td></tr><tr><td>Effective temperature, T<sub>eff</sub></td><td>K</td><td>6407<sup>+81</sup><sub>-78</sub></td></tr><tr><td>Radius, R<sub>s</sub></td><td>R<sub>☉</sub></td><td>1.49 ± 0.25</td></tr><tr><td>Metallicity, [M/H]</td><td>dex</td><td>-0.02 ± 0.10</td></tr><tr><td>Extinction, A<sub>V,S</sub></td><td></td><td>1.11 ± 0.05</td></tr><tr><td>Angular radius, θ<sub>*</sub></td><td>μas</td><td>8.63 ± 0.06</td></tr><tr><td>Distance, D<sub>s</sub></td><td>10<sup>2</sup> pc</td><td>8.0 ± 1.3</td></tr></tbody></table>

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 ($\theta_*$) 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_{\text{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_{\text{eff}}$, $R_S$, [M/H], and $A_{V,S}$ by adding penalties to the $\chi^2$ value as

$$ \chi^2 = \sum_{\lambda} \frac{(f_{\text{obs},\lambda} - f_{\text{model},\lambda})^2}{\sigma_{f_{\text{obs},\lambda}}^2} + \sum_{i} \frac{(X_i - X_{i,\text{prior}})^2}{\sigma_{X_{i,\text{prior}}}^2}, \quad (7) $$

where $f_{\text{obs},\lambda}$, $\sigma_{f_{\text{obs},\lambda}}$, and $f_{\text{model},\lambda}$ are the observed flux density, its $1\sigma$ uncertainty, and the model flux density, respectively, for a band $\lambda$, and $X_i$ denotes one of the parameters among $T_{\text{eff}}$, $R_S$, [M/H], and $A_V$. For the priors of $T_{\text{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\sigma$ uncertainty.

The derived median value and $1\sigma$ 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 \pm 130$ pc and $\theta_* = 8.63 \pm 0.06$ µas, respectively, which are well consistent with the previous estimations of $D_S = 700-800$ pc (Nucita et al. 2018) and $\theta_* = 9 \pm 0.9$ µas (Dong et al. 2019).

# 5. Physical Parameters of the Lens System

## 5.1. Constraint from the Microlensing Model

If $\theta_E$, $\pi_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:

$$ M_L = \frac{\theta_E}{\kappa\pi_E}, \quad (8) $$

$$ D_L = \frac{AU}{\pi_E \theta_E + \pi_S}, \quad (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_{\text{proj}} = s\theta_E D_L$. The median and $1\sigma$ uncertainties of these parameters derived from the light-curve analysis using the $\theta_E$ and $\Phi_\pi$ 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.

However, as discussed in Section 3.5.1, the detection of $\pi_E$ is marginal, and the signal is as weak as the level of systematics. Therefore, it is conservative not to rely on the $\pi_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. From the Lens Brightness

### 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.

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'' \times 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 $\theta_E$ prior and the upper limit on $\pi_E$ (<1.1), we calculate the probability distributions of the fraction of the lens flux to the total blending flux, $f_L \equiv 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 $\theta_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_\odot$, 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
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Corner plot for the parameters of the source star. The black and gray areas indicate the 68% and 95% confidence regions, respectively. Note that the bimodal feature in [M/H] centered at [M/H] = 0 is an artifact due to the discreteness of the theoretical models we adopt.</figcaption></figure>

$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

$$R_{*} = \sum_{i}^{n} a_{i} M_{\lambda}^{i} \times (1 + f [\text{Fe/H}]), \quad (10)$$

where $R_*$ is the stellar radius, $M_\lambda$ is the absolute magnitude in the $\lambda$ 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

$$m_{\lambda, \text{calc}} = M_{\lambda} + 5 \log_{10}(D_{L}/10\text{pc}) + A_{\lambda, L}, \quad (11)$$

where $\lambda$ 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 $\lambda$ band. Note that $M_\lambda$ is tied with the radius, $R_{L1}$, and metallicity, [Fe/H], of the lens star via Equation (10). Here we adopt

$A_\lambda/A_V = (1.223, 1.011, 0.880, 0.676, 0.485, 0.117)$ for $\lambda = (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 $\chi^2$ value:

$$\chi^2 = \sum_{\lambda=\{g', r', V, i', z', K_s\}} \frac{(m_{\lambda, \text{obs}} - m_{\lambda, \text{calc}})^2}{\sigma^2_{m_{\lambda, \text{obs}}}} + \frac{([\text{Fe/H}] - [\text{Fe/H}]_{\text{prior}})^2}{\sigma^2_{[\text{Fe/H}]_{\text{prior}}}} + \frac{(A_{V,L} - A_{V,L, \text{prior}})^2}{\sigma^2_{A_{V,L, \text{prior}}}}, \quad (12)$$

where $m_{\lambda, \text{obs}}$ and $\sigma_{m_{\lambda, \text{obs}}}$ are the observed magnitude and its $1\sigma$ uncertainty in the $\lambda$ band, respectively; $[\text{Fe/H}]_{\text{prior}}$ is a prior for [Fe/H]; and $A_{V,L, \text{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, \text{prior}}$ and 0.10 to $\sigma_{A_{V,L, \text{prior}}}$ in the same way as for $A_{V,s}$. The derived
<table><caption>**Table 4**

Physical Parameters of the Lens System</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Nucita et al. (2018)</th><th>π<sub>E</sub> and θ<sub>E</sub> and D<sub>S</sub></th><th>Lens Flux</th><th>Lens Flux and θ<sub>E</sub> and D<sub>S</sub></th></tr></thead><tbody><tr><td>Distance, D<sub>L</sub></td><td>pc</td><td>~380</td><td>511<sup>+101</sup><sub>-80</sub></td><td>507 ± 74</td><td>505 ± 47</td></tr><tr><td>Stellar mass, M<sub>L1</sub></td><td>M<sub>☉</sub></td><td>0.25 ± 0.18</td><td>0.64<sup>+0.38</sup><sub>-0.19</sub></td><td>0.590<sup>+0.042</sup><sub>-0.051</sub></td><td>0.586 ± 0.033</td></tr><tr><td>Stellar radius, R<sub>L1</sub></td><td>R<sub>☉</sub></td><td>...</td><td>...</td><td>0.599<sup>+0.056</sup><sub>-0.061</sub></td><td>...</td></tr><tr><td>Extinction, A<sub>V,L</sub></td><td></td><td>...</td><td>...</td><td>0.95 ± 0.11</td><td>...</td></tr><tr><td>Metallicity, [Fe/H]</td><td>dex</td><td>...</td><td>...</td><td>-0.05 ± 0.20</td><td>...</td></tr><tr><td>Absolute K<sub>s</sub> magnitude, M<sub>Ks</sub></td><td>mag</td><td>...</td><td>...</td><td>5.05<sup>+0.33</sup><sub>-0.28</sub></td><td>...</td></tr><tr><td>Planetary mass, M<sub>L2</sub></td><td>M<sub>☉</sub></td><td>9.2 ± 6.6</td><td>21.8<sup>+12.9</sup><sub>-6.5</sub></td><td>20.0 ± 2.3</td><td>20.0 ± 2.0</td></tr><tr><td>Projected separation, a<sub>proj</sub> (model a)</td><td>au</td><td>~0.5</td><td>0.89<sup>+0.18</sup><sub>-0.14</sub></td><td>0.89 ± 0.13</td><td>0.88 ± 0.08</td></tr><tr><td>Projected separation, a<sub>proj</sub> (model b)</td><td>au</td><td>~0.5</td><td>0.95<sup>+0.19</sup><sub>-0.15</sub></td><td>0.95 ± 0.14</td><td>0.94 ± 0.09</td></tr><tr><td>Semimajor axis, a<sub>circ</sub><sup>a</sup></td><td>au</td><td>...</td><td>1.12<sup>+0.66</sup><sub>-0.25</sub></td><td>1.10<sup>+0.63</sup><sub>-0.22</sub></td><td>1.08<sup>+0.62</sup><sub>-0.18</sub></td></tr></tbody></table>

**Note.**

<sup>a</sup> Calculated by merging the posteriors of models *a* and *b*.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Posterior distributions of the mass and distance of the lens star. Blue dotted contours, gray shaded regions, red solid contours, and green shaded regions indicate the constraints calculated from $\pi_E$, $\theta_E$, and $D_S$; $\theta_E$ and $D_S$; lens flux; and the combination of lens flux and $\theta_E$ and $D_S$, respectively. In each case, dark (inner) and light (outer) colored lines or shaded regions represent 68% and 95% confidence regions, respectively. The cyan dashed line indicates a lower limit given by the 3$\sigma$ upper limit of $\pi_E$ and 3$\sigma$ lower limit of $D_S$.</figcaption></figure>

<table><caption>**Table 5**

Calibrated Magnitudes of the Blending Flux</caption><thead><tr><th>Band</th><th>Magnitude</th></tr></thead><tbody><tr><td>g'</td><td>19.088 ± 0.337</td></tr><tr><td>V</td><td>17.760 ± 0.110</td></tr><tr><td>r'</td><td>17.305 ± 0.122</td></tr><tr><td>i'</td><td>16.382 ± 0.068</td></tr><tr><td>z'</td><td>15.872 ± 0.051</td></tr><tr><td>K<sub>s</sub></td><td>13.728 ± 0.027</td></tr></tbody></table>

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\sigma$ 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.

### 5.3. Combined Solution

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_E$ and $D_S$ instead of the posterior distribution of the $M_{L1}$ and $D_L$ solution from $\pi_E$, $\theta_E$, and $D_S$, because the latter relies on the posterior distribution of $\pi_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_{proj} \equiv s\theta_E D_L = 0.88 \pm 0.08$ au (model *a*) and $0.94 \pm 0.09$ au (model *b*), which is converted to the semimajor axis of $a_{circ} = 1.08^{+0.62}_{-0.18}$ au, where a circular orbit and random orientation are assumed and the solutions of two models (model *a* and *b*) are merged.

## 6. Discussions

### 6.1. Comparison of the Planetary Location with the Snow Line

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),
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Corner plot for the parameters of the lens star derived from the lens brightness. The black and gray areas indicate the 68% and 95% confidence regions, respectively.</figcaption></figure>

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%.

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_{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.

According to the core accretion theories, it is difficult to form a planet as massive as Kojima-1Lb ($20 \pm 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; Drąż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
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Distribution of known exoplanets in the planetary mass and semimajor axis planes for the host stars having a mass of 0.4–0.8 $M_\odot$. Data are collected mainly from http://exoplanet.eu. Black squares, blue circles, and red circles indicate the planets observed by radial velocity, transit, and microlensing, respectively. The filled and open circles of microlensing show the planets with and without direct mass constraint, respectively. Two degenerated solutions are connected by a dotted line, if applicable. Kojima-1Lb is depicted as a green circle. The contours show the planet detection efficiencies for Kojima-1 of 90%, 70%, 40%, and 10% (top to bottom). (b) Same as panel (a), but the x-axis is converted to the semimajor axis normalized by the snow-line location estimated by $a_{\text{snow}} = 2.7 \times M_\odot$ au.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Snow-line distance as a function of time. The solid line indicates a theoretical model for a disk of a $0.6 M_\odot$ star considering stellar irradiation and viscous accretion, extracted from Figure 1 of Kennedy & Kenyon (2008). The dashed line is a time-independent snow-line location for Kojima-1L calculated by $a_{\text{snow}} = 2.7 \times M_\odot/M_\odot$ au. The median value and $1\sigma$ confidence region of the semimajor axis of Kojima-1Lb are shown as a gray dotted line and light gray shaded area, respectively.</figcaption></figure>

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.

## 6.2. Detection Efficiency to the Planetary Signal

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 $\Delta\chi^2 = 100$ following Suzuki et al. (2016), where $\Delta\chi^2$ is the $\chi^2$ difference between planetary and nonplanetary (single-lens) models. At first, the detection efficiency $\epsilon$ is computed as a function of (log s, log q). Next, we transform it to the physical parameter space, (log $a_{\text{proj}}$, 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$ (log $a_{\text{proj}}$, log $M_{L2}$) is further converted to $\epsilon$ (log $a_{3D}$, 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 log($a_{3D}/a_{\text{snow}}$) and log $M_{L2}$, where $a_{\text{snow}} = 2.7 \times (M_\odot/M_\odot)$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.
## 6.3. Capabilities of Future Follow-up Observations

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_{geo} = 25.34 \pm 0.44$ mas yr$^{-1}$, enabling us to spatially separate the source and lens stars in $\sim 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 $\sim 2.2$ sin $i$ ms$^{-1}$ and a period of $\sim 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.

## 7. Summary

We conducted follow-up observations of the nearby planetary microlensing event Kojima-1 by means of seeing-limited photometry, spectroscopy, and high-resolution imaging. We found no additional planetary feature in our photometric data other than the one that was identified by Nucita et al. (2017). From the light-curve modeling and spectroscopic analysis, we have refined the distance and angular diameter of the source star to be 800 $\pm$ 130 pc and 8.63 $\pm$ 0.06 $\mu$as, respectively. We have also refined the microlensing model using the prior information of $\theta_E$ and $\Phi_\pi$ from the VLTI observation by Dong et al. (2019). We confirm the presence of apparent blending flux and absence of significant parallax signal reported in the literature. We find no contaminating sources in the Keck AO image and that the detected blending flux most likely comes from the lens star. Combining all of this information, we have directly derived the physical parameters of the lens system without relying on any Galactic models, finding that the host star is a dwarf on the M/K boundary ($0.59 \pm 0.03 M_{\odot}$) located at 500 $\pm$ 50 pc, and the companion is a Neptune-mass planet ($20 \pm 2 M_{\oplus}$) with a semimajor axis of $\sim 1.1$ au.

The orbit of Kojima-1Lb is a few times closer to the host star than the other microlensing planets around the same type of star

and is likely comparable to the snow-line distance at its youth. We have estimated that the detection efficiency of this planet in this event is $\sim 35\%$, which may imply that Neptunes are common around the snow line.

The host star is the brightest ($K_s = 13.7$) among all of the microlensing planetary systems, providing us a great opportunity not only to spectroscopically characterize the host star but also to confirm the mass and refine the orbit of this planet with the radial velocity technique in the near future.

We thank the anonymous referee for a lot of thoughtful comments. A.F. thanks T. Kimura and H. Kawahara for meaningful discussions on the formation and abundance of Neptunes around the snow line. A.F. also thanks A. Nucita and A. Mann for kindly providing data used in their papers.

This article is based on observations made with the MuSCAT2 instrument, developed by ABC, at Telescopio Carlos Sánchez, operated on the island of Tenerife by the IAC in the Spanish Observatorio del Teide. We acknowledge ISAS/JAXA for the use of its facility through the inter-university research system. A.Y. is grateful to Mizuki Isogai, Akira Arai, and Hideyo Kawakita for their technical support on observations with the Araki telescope. D.S. acknowledges The Open University for the use of the COAST telescope. A.F. acknowledges the MOA collaboration/Osaka University for the use of the computing cluster.

This work was partly supported by JSPS KAKENHI grant Nos. JP25870893, JP16K17660, JP17H02871, JP17H04574, JP18H01265, and JP18H05439; MEXT KAKENHI grant Nos. JP17H06362 and JP23103004; and JST PRESTO grant No. JPMJPR1775. This work was also partially supported by the Optical and Near-Infrared Astronomy Inter-University Cooperation Program of the MEXT of Japan and the JSPS and NSF under the JSPS-NSF Partnerships for International Research and Education. This work was partly financed by the Spanish Ministry of Economics and Competitiveness through grants ESP2013-48391-C4-2-R and AYA2015-69350-C3-2-P. Y.T. acknowledges the support of DFG priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (WA 1047/11-1). S.Sh. acknowledges support from grants APVV-15-0458 and VEGA 2/0008/17.

## Appendix

To complement Table 1 of Mann et al. (2015), we calculate the coefficients of the radius–metallicity–luminosity relation for other bands than $K_s$ using the same data set used by Mann et al. (2015). They made public a table that includes synthetic apparent magnitudes in various bands (calculated from cataloged magnitudes and low-resolution spectra) and stellar radius (estimated from the observed bolometric flux and effective temperature) for 183 nearby M7–K7 single stars. This table, however, lacks the information on parallax that is needed to convert the apparent magnitude to absolute magnitude, which we got from the authors by private communication. (Their parallax came from somewhere before Gaia, but we do not attempt to update them using Gaia to keep consistency.)

To derive the relation, we apply Equation (5) of their paper, that is,

$$R_{*} = (a + bM_{\lambda} + cM_{\lambda}^{2} + \dots) \times (1 + f[\text{Fe/H}]), \quad (13)$$

where $R_{*}$ is the stellar radius, $M_{\lambda}$ is the absolute magnitude in band $\lambda$, [Fe/H] is the metallicity, and a, b, c, ..., f are
<table><caption>**Table 6**
Coefficients of Radius–Metallicity–Luminosity Relation</caption><thead><tr><th>Band</th><th>a</th><th>b</th><th>c</th><th>d</th><th>e</th><th>f</th></tr></thead><tbody><tr><td>g'</td><td>-4.0294</td><td>1.6103</td><td>-1.9349 × 10<sup>-1</sup></td><td>9.4899 × 10<sup>-3</sup></td><td>-1.6655 × 10<sup>-4</sup></td><td>3.2209 × 10<sup>-1</sup></td></tr><tr><td>r'</td><td>-2.5349</td><td>1.2698</td><td>-1.7485 × 10<sup>-1</sup></td><td>9.6309 × 10<sup>-3</sup></td><td>-1.8821 × 10<sup>-4</sup></td><td>3.4127 × 10<sup>-1</sup></td></tr><tr><td>i'</td><td>-3.5485</td><td>1.9081</td><td>-2.9955 × 10<sup>-1</sup></td><td>1.9070 × 10<sup>-2</sup></td><td>-4.3370 × 10<sup>-4</sup></td><td>2.5015 × 10<sup>-1</sup></td></tr><tr><td>z'</td><td>-3.9416</td><td>2.3156</td><td>-4.0010 × 10<sup>-1</sup></td><td>2.8101 × 10<sup>-2</sup></td><td>-7.0665 × 10<sup>-4</sup></td><td>1.766 × 10<sup>-1</sup></td></tr><tr><td>V</td><td>-3.1842</td><td>1.4307</td><td>-1.8538 × 10<sup>-1</sup></td><td>9.7067 × 10<sup>-3</sup></td><td>-1.8107 × 10<sup>-4</sup></td><td>3.3462 × 10<sup>-1</sup></td></tr><tr><td>K<sub>s</sub></td><td>1.9305</td><td>-3.4665 × 10<sup>-1</sup></td><td>1.6472 × 10<sup>-2</sup></td><td>...</td><td>...</td><td>4.4889 × 10<sup>-2</sup><sup>a</sup></td></tr></tbody></table>

**Note.**

ª There is a small difference in the values between this work and Mann et al. (2015), which we suspect due to round errors in [Fe/H].

coefficients. We choose the polynomial order for $M_\lambda$ such that the best-fit BIC value (Schwarz 1978) is minimized. We derive the coefficients for the $g'$, $r'$, $i'$, $z'$, and $V$ bands, as well as for the $K_s$ band, for completeness, as listed in Table 6.

**ORCID iDs**

A. Fukui https://orcid.org/0000-0002-4909-5763

E. Bachelet https://orcid.org/0000-0002-6578-5078

D. Mawet https://orcid.org/0000-0002-8895-4735

M. T. Penny https://orcid.org/0000-0001-7506-5640

N. Narita https://orcid.org/0000-0001-8511-2981

M. Fukagawa https://orcid.org/0000-0003-3500-2455

Y. Tsapras https://orcid.org/0000-0001-8411-351X

M. Hundertmark https://orcid.org/0000-0003-0961-5231

D. E. Reichart https://orcid.org/0000-0002-5060-3673

M. Tamura https://orcid.org/0000-0002-6510-0681

**References**

Alard, C. 2000, A&AS, 144, 363

Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017, A&A, 602, A88

Bachelet, E., Norbury, M., Bozza, V., & Street, R. 2017, AJ, 154, 203

Batista, V., Beaulieu, J.-P., Bennett, D. P., et al. 2015, ApJ, 808, 170

Bennett, D. P. 2010, ApJ, 716, 1408

Bennett, D. P., Bhattacharya, A., Anderson, J., et al. 2015, ApJ, 808, 169

Bennett, D. P., Bond, I. A., Udalski, A., et al. 2008, ApJ, 684, 663

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bhattacharya, A., Beaulieu, J.-P., Bennett, D. P., et al. 2018, AJ, 156, 289

Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112

Bozza, V., Bachelet, E., Bartolić, F., et al. 2018, MNRAS, 479, 5157

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77

Dominik, M. 2006, MNRAS, 367, 669

Dong, S., Mérand, A., Delplancke-Ströbele, F., et al. 2019, ApJ, 871, 70

Drążkowska, J., & Alibert, Y. 2017, A&A, 608, A92

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fukui, A., Abe, F., Ayani, K., et al. 2007, ApJ, 670, 423

Fukui, A., Narita, N., Tristram, P. J., et al. 2011, PASJ, 63, 287

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Gaidos, E., & Mann, A. W. 2014, ApJ, 791, 54

Gaudi, B. S., Patterson, J., Spiegel, D. S., et al. 2008, ApJ, 677, 1268

Han, C. 2008, ApJ, 681, 806

Han, C., Yee, J. C., Udalski, A., et al. 2019, AJ, 158, 102

Hansen, B. M. S., & Murray, N. 2012, ApJ, 751, 158

Herman, M. K., Zhu, W., & Wu, Y. 2019, AJ, 157, 248

Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388

Ida, S., & Lin, D. N. C. 2005, ApJ, 626, 1045

Jayasinghe, T., Dong, S., Stanek, K. Z., et al. 2017, ATel, 10923, 1

Kambe, E., Yoshida, M., Izumiura, H., et al. 2013, PASJ, 65, 15

Kawahara, H., & Masuda, K. 2019, AJ, 157, 218

Kennedy, G. M., & Kenyon, S. J. 2008, ApJ, 673, 502

Kennedy, G. M., Kenyon, S. J., & Bromley, B. C. 2006, ApJL, 650, L139

Kesseli, A. Y., West, A. A., Veyette, M., et al. 2017, ApJS, 230, 16

Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502

Kokubo, E., & Ida, S. 2002, ApJ, 581, 666

Konyves-Toth, R., Pal, A., Ordasi, A., & Vinko, J. 2017, ATel, 10926, 1

Koshimoto, N., Shvartzvald, Y., Bennett, D. P., et al. 2017, AJ, 154, 3

Kotani, T., Kawai, N., Yanagisawa, K., et al. 2005, NCIMC, 28, 755

Luhman, K. L., & Mamajek, E. E. 2012, ApJ, 758, 31

Maehara, H. 2017, ATel, 10919, 1

Mann, A. W., Dupuy, T., Kraus, A. L., et al. 2019, ApJ, 871, 63

Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64

Mordasini, C., Alibert, Y., Benz, W., & Naef, D. 2009, A&A, 501, 1161

Mordasini, C., Mayor, M., Udry, S., et al. 2011, A&A, 526, A111

Muraki, Y., Han, C., Bennett, D. P., et al. 2011, ApJ, 741, 22

Narita, N., Fukui, A., Kusakabe, N., et al. 2015, JATIS, 1, 045001

Narita, N., Fukui, A., Kusakabe, N., et al. 2019, JATIS, 5, 015001

Nucita, A. A., Licchelli, D., De Paolis, F., et al. 2018, MNRAS, 476, 2962

Nucita, A. A., Licchelli, D., De Paolis, F., Ingrosso, G., & Strafella, F. 2017, ATel, 10934, 1

Ogihara, M., Morbidelli, A., & Guillot, T. 2015, A&A, 578, A36

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

Rhie, S. H., Bennett, D. P., Becker, A. C., et al. 2000, ApJ, 533, 378

Ribas, Á., Bouy, H., & Merín, B. 2015, A&A, 576, A52

Schwarz, G. 1978, AnSta, 6, 461

Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48

Sokolovsky, K. 2017, ATel, 10921, 1

Sumi, T., Bennett, D. P., Bond, I. A., et al. 2010, ApJ, 710, 1641

Suzuki, D., Bennett, D. P., Ida, S., et al. 2018, ApJL, 869, L34

Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145

Torres, G., Andersen, J., & Giménez, A. 2010, A&ARv, 18, 67

Tsapras, Y. 2018, Geosc, 8, 365

Tuomi, M., Jones, H. R. A., Butler, R. P., et al. 2019, arXiv:1906.04644

Udalski, A., Ryu, Y.-H., Sajadian, S., et al. 2018, AcA, 68, 1

Villanueva, S., Jr., Gaudi, B. S., Pogge, R. W., et al. 2018, PASP, 130, 015001

Yanagisawa, K., Kuroda, D., Yoshida, M., et al. 2010, in AIP Conf. Ser. 1279, Deciphering the Ancient Universe with Gamma-Ray Bursts, ed. N. Kawai & S. Nagataki (Melville, NY: AIP), 466

Yanagisawa, K., Shimizu, Y., Okita, K., et al. 2016, Proc. SPIE, 9908, 99085D

Yee, S. W., Petigura, E. A., & von Braun, K. 2017, ApJ, 836, 77