Four sub-Jovian-mass planets detected by high-cadence microlensing surveys

With the aim of finding short-term planetary signals, we investigated the data collected from the high-cadence microlensing surveys. From this investigation, we found four planetary systems with low planet-to-host mass ratios, including OGLE-2017-BLG-1691L, KMT-2021-BLG-0320L, KMT-2021-BLG-1303L, and KMT-2021-BLG-1554L. Despite the short durations, ranging from a few hours to a couple of days, the planetary signals were clearly detected by the combined data of the lensing surveys. It is found that three of the planetary systems have mass ratios of the order of $10^{-4}$ and the other has a mass ratio slightly greater than $10^{-3}$. The estimated masses indicate that all discovered planets have sub-Jovian masses. The planet masses of KMT-2021-BLG-0320Lb, KMT-2021-BLG-1303Lb, and KMT-2021-BLG-1554Lb correspond to $\sim 0.10$, $\sim 0.38$, and $\sim 0.12$ times of the mass of the Jupiter, and the mass of OGLE-2017-BLG-1691Lb corresponds to that of the Uranus. The estimated mass of the planet host KMT-2021-BLG-1554L, $M_{\rm host}\sim 0.08~M_\odot$, corresponds to the boundary between a star and a brown dwarf. Besides this system, the host stars of the other planetary systems are low-mass stars with masses in the range of $\sim [0.3$--$0.6]~M_\odot$. The discoveries of the planets well demonstrate the capability of the current high-cadence microlensing surveys in detecting low-mass planets.


Introduction
During the 2010s, microlensing entered an era of high-cadence surveys with the instrumental upgrade of previously established experiments and the participation of a new experiment. The new era started when the Microlensing Observations in Astrophysics survey (MOA: Bond et al. 2001), which had previously carried out a lensing survey using a 0.61 m telescope and a camera with a 1.3 deg 2 field of view (FOV) during the early phase, launched its second phase experiment with the employment of a 1.8 m telescope equipped with a wide-field camera yielding a 2.2 deg 2 FOV. Since its first operation in 1992, the Optical Gravitational Lensing Experiment (OGLE) has been upgraded multiple times, and it is in its fourth phase (OGLE-IV: Udalski et al. 2015) with the employment of a 1.3 m telescope and a camera with a 1.4 deg 2 FOV. The Korea Microlensing Telescope Net-work (KMTNet: Kim et al. 2016), which was launched in 2016, is being carried out with the use of three globally distributed 1.6 m telescopes, each of which is mounted with a wide-field detector providing a 4 deg 2 FOV. With the use of wide-field cameras mounted on multiple telescopes, the current lensing surveys can monitor lensing events with a dramatically enhanced observational cadence.
The greatly increased observational cadence of the lensing surveys has brought out important changes not only in the observational strategy but also in the outcome of planetary microlensing searches. First, being able to resolve short-lasting planetary signals from survey observations, the high-cadence surveys can function without the survey+followup experiment mode, in which low-cadence surveys mainly detect and alert lensing events and followup groups densely observe the alerted events. Second, the number of microlensing planets has substan- Notes. HJD ′ = HJD − 2450000. tially increased with the operation of the high-cadence surveys. Among the total 135 microlensing planets with measured mass ratios 1 , 80 (59%) planets were found with the use of the data from the KMTNet survey. This can be seen in Figure 1, in which we present the histogram of microlensing planets as a function of the planet-to-host mass ratio q. The histogram also shows that the high-cadence surveys contribute to the detections of planets with very low mass ratios, especially those with q < 10 −4 . In this paper, we report four planetary systems with low planet-to-host mass ratios detected from the microlensing surveys. Despite the short durations, ranging from a few hours to a 1 The sample includes 119 planets listed in the NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/index.html) plus 16 planets that are not included in the archive: KMT-2020-BLG-0414Lb (Zang et al. 2021a), OGLE-2019-BLG-1053Lb (Zang et al. 2021b), KMT-2019-BLG-0253Lb, KMT-2019-BLG-0953Lb, OGLE-2018-BLG-0977Lb, OGLE-2018-BLG-0506Lb, OGLE-2018-BLG-0516Lb, OGLE-2019-BLG-1492Lb , KMT-2017-BLG-2509Lb, OGLE-2017-BLG-1099Lb, OGLE-2019-BLG-0299Lb ), OGLE-2016-BLG-1093Lb (Shin et al. 2022, KMT-2018-BLG-1988Lb (Han et al. 2022b), KMT-2021-BLG-0912Lb (Han et al. 2022a), OGLE-2019-BLG-0468Lb, and OGLE-2019-BLG-0468Lc (Han et al. 2022c). couple of days, the planetary signals were clearly detected by the combined data of the lensing surveys. It is found that all of the discovered planets have sub-Jovian masses lying in the range of [0.05-0.38] M J , illustrating the importance of the high-cadence surveys in detecting low-mass planets.
We present the analysis of the planetary lensing events according to the following organization. In Sect. 2, we describe the observations of the individual lensing events, instrument used for observations, and the process of the data reduction. In Sect. 3, we explain the procedure of analysis that is applied to each of the individual events. In the subsequent subsections, we describe the detailed features of the planetary signals and present the results of the analyses conducted for the individual lensing events. In Sect. 4, we specify the types of the source stars and estimate the angular Einstein radii of the lens systems. In Sect. 5, we estimate the physical parameters of the planetary systems by conducting Bayesian analyses. We summarize the results of the analyses and conclude in Sect. 6.
The events were detected from the combined observations of the lensing surveys conducted by the KMTNet, OGLE, and MOA groups. The KMTNet survey utilizes three identical 1.6 m telescopes that are distributed in three countries of the Southern Hemisphere for the continuous coverage of lensing events. The sites of the individual telescopes are the Siding Spring Observatory in Australia (KMTA), the Cerro Tololo InterAmerican Observatory in Chile (KMTC), and the South African Astronomical Observatory in South Africa (KMTS). Observations by the OGLE survey were conducted using the 1.3 m telescope located at Las Campanas Observatory in Chile. The MOA survey uses the 1.8 m telescope at the Mt. John Observatory in New Zealand. Observations were done mainly in the I band for the KMTNet and OGLE surveys and in the customized MOA-R band for the MOA survey. A subset of images were acquired in the V band for the purpose of estimating source colors. The detailed procedure of the source color estimation will be discussed in Sect. 4.
Reduction and photometry of the data were done using the pipelines of the individual survey groups: Albrow et al. (2009) for the KMTNet survey, Woźniak (2000) for the OGLE survey, and Bond et al. (2001) for the MOA survey. All these pipelines apply the difference image analysis algorithm  (Albrow 2017) for the specification of the source colors. Following the routine described in Yee et al. (2012), we rescale the error bars of data by σ = k(σ 2 min + σ 2 0 ) 1/2 , where σ 0 represents the error estimated from the photometry pipeline, σ min is a factor used to make the data consistent with the scatter of data, and the factor k is used to make χ 2 per degree of freedom for each data set become unity. In Table 2, we list the factors k and σ min of the individual data sets.

Analyses
The light curves of the analyzed events share a common characteristic that short-lived anomalies appear on the otherwise smooth and symmetric form of a single-lens single-source (1L1S) event. Such anomalies in lensing light curves can be produced by two channels, in which the first channel is a perturbation induced by a low-mass companion, such as a planet, to the lens (Gould & Loeb 1992), and the second channel is an anomaly caused by a faint companion to the source (Gaudi 1998). Hereafter we denote events with a binary lens and a binary source as 2L1S and 1L2S events, respectively. We examine the origins of the anomalies by modeling the light curve of the individual events under these 2L1S and 1L2S interpretations.
Lensing light curves are described by the combination of various lensing parameters. For a 1L1S event, the light curve is described by three parameters of (t 0 , u 0 , t E ), which denote the time of the closest approach of the source to the lens, the separation (scaled to the angular Einstein radius θ E ) between the source and lens at that time, and the time scale of the event, respectively. The event time scale is defined as the time required for a source to transit θ E .
In addition to these basic parameters, describing the light curve of a 2L1S event requires one to include extra parameters of (s, q, α), which represent the separation (scaled to θ E ) and mass ratio between the binary lens components, and the angle between the binary-lens axis and the direction of the source motion (source trajectory angle), respectively. Because planetinduced anomalies are usually generated by the crossings over or close approach of the source to the planet-induced caustics, an additional parameter of the normalized source radius ρ, which is defined as the ratio of the angular radius of the source θ * to θ E , is needed to describe the deformation of the anomaly by finite-source effects (Bennett & Rhie 1996). In computing finitesource magnifications, we consider the limb-darkening variation of the source.
One also needs extra parameters to describe the lensing light curve of a 1L2S event. These parameters are (t 0,2 , u 0,2 , q F ), which represent the closest approach time and separation between the lens and the second source, and the flux ratio between the binary source stars, respectively. In the 1L2S model, we designate the lensing parameters related to the primary source, S 1 , as (t 0,1 , u 0,1 ) to distinguish them from those related to the source companion, S 2 . In order to consider finite-source effects occurring when the lens passes over either of the source stars, we add two additional parameters of the normalized source radii ρ 1 and ρ 2 for the lens transits over S 1 and S 2 , respectively (Dominik 2019).
In modeling 1L1S and 1L2S light curves, for which the lensing magnification is smooth with the variation of the lensing parameters, we search for the solution of the lensing parameters via a downhill approach by minimizing χ 2 using the Markow Chain Monte Carlo (MCMC) algorithm. In modeling 2L1S light curves, for which the lensing magnification variation is discontinuous due to the formation of caustics and there may exist multiple local solutions caused by various types of degeneracy, we model the light curves in two steps. In the first step, we conduct grid searches for the binary parameters s and q, construct a χ 2 map on the s-q parameter plane, and then identify local solutions appearing on the χ 2 map. In the second step, we polish the individual local solutions using the downhill approach. Below we present the details of the modeling conducted for the individual lensing events.

OGLE-2017-BLG-1691 (KMT-2017-BLG-0752)
The lensing event OGLE-2017-BLG-1691 occurred during the 2017 bulge season. It was first found by the OGLE survey, and later confirmed by the KMTNet survey from the post-season analyses of the data collected during the season. The KMTNet group designated the event as KMT-2017-BLG-0752. Hereafter, we use the nomenclatures of events by the ID references of the surveys who first found the events in accordance with the convention of the microlensing community. The baseline magnitude of the event was I base = 19.900 ± 0.004.
The lensing light curve of OGLE-2017-BLG-1691 is shown in Figure 2, in which the bottom panel shows the whole view and the top panel displays the enlargement of the region around the anomaly. The event was alerted by the OGLE survey on 2017 September 6, HJD ′ ≡ HJD − 2450000 ∼ 8002.5, which approximately corresponds to the time near the peak. An anomaly occurred about one day after the peak, but it was not noticed during the progress of the event because it was covered by just a single OGLE data point, and the KMTNet data were not released during the lensing magnification. The existence of the anomaly was identified five years after the event, shortly after the re-reduction of all 2017 KMT light curves, from a project of reinvestigating the previous KMTNet data conducted to find unnoticed planetary signals (Han et al. 2022b). The top panel of Figure 2 shows that the anomaly, which lasted for ∼ 5 hours centered at HJD ′ ∼ 8003.6, was additionally covered by three KMTC data points, and this confirmed that the single anomalous OGLE point was real.
Notes. HJD ′ = HJD − 2450000. In the top panel, curves of the three tested models, 2L1S (outer), 1L2S, and 1L1S models, are drawn over the data points, and the residuals from the individual models are presented in the three middle panels. The curve drawn in the bottom panel is the best-fit model (outer 2L1S model). Colors of data points are set to match those of the telescopes used for observations marked in the legend. The curves drawn in the 1L2S and 1L1S residual panels represent the differences from the 2L1S model.
We model the light curve of the event under the 1L1S, 2L1S and 1L2S interpretations. The residuals of the individual tested models are compared in Figure 2. It is found that the 2L1S model best describes the anomaly, being favored over the 1L1S and 1L2S models by ∆χ 2 = 456.1 and 13.9, respectively. From the comparison of the residuals, it is found that the 2L1S model is confirmed not only by the 4 data points (3 KMTC plus 1 OGLE points) with positive deviations near the peak of the bump but also by the 3 additional points (2 KMTS and 1 KMTC points) before the major bump with slightly negative deviations. We find two degenerate 2L1S solutions with binary parameters of (s, q) ∼ (1.058, 0.71 × 10 −4 ) and ∼ (1.003, 0.97 × 10 −4 ), which we designate as "inner" and "outer" solutions, respectively, for the reason to be mentioned below. We list the full lensing parameters of the two 2L1S solutions in Table 3, together with the parameters of the 1L2S model. The degeneracy is very severe and the outer solution is favored only by ∆χ 2 = 0.4. The lensing parameters of the solutions indicate that the anomaly was produced by a very small mass-ratio planetary companion to the lens lying close to the Einstein ring of the primary regardless of the solutions. From the fact that both solutions have separations greater than unity and do not follow the relation of s inner × s outer ≃ 1, the degeneracy is different from the "close-wide" degeneracy (Griest & Safizadeh 1998;Dominik 1999;An 2005), which arises due to the similarity between the central caustics of planetary lens systems with planet separations s and 1/s. Instead, the planet separations of the two degenerate solutions follow the relation where u 2 anom = τ 2 anom + u 2 0 , τ anom = (t anom − t 0 )/t E , and t anom indicates the time of the planetary anomaly Zhang et al. 2022;Ryu et al. 2022). With the values t 0 ≃ 8002.5, u 0 ≃ 4.85 × 10 −2 , t E ≃ 19 day and t anom ≃ 8003.6, one finds that s † ≃ 1.04, which matches well √ s in × s out ≃ 1.03. This indicates that the similarity between the light curves of the two solutions is caused by the degeneracy identified by Yee et al. (2021), who first mentioned the continuous transition between the "close-wide" and "inner-outer" (Gaudi & Gould 1997) degeneracies. Hereafter, we refer to this degeneracy as "offset degeneracy" following Zhang et al. (2022). Figure 3 shows the lens systems configurations, in which the source trajectory (line with an arrow) with respect to the caustic and positions of the lens components (blue solid dots marked by M 1 and M 2 ) is presented: inner solution in the upper panel and outer solution in the lower panel. The source passed the inner side (with respect to M 1 ) of the caustic according to the inner solution, while the source passes the outer side according to the outer solution. For both solutions, the source crossed the cusp of the caustic, and thus the anomaly is affected by finitesource effects during the caustic crossing, allowing us to precisely measure the normalized source radius. For each solution, the source size relative to the caustic is represented by an orange circle marked on the source trajectory. As will be discussed in Sect. 4, the measurement of ρ is important to estimate the lensing observable of the angular Einstein radius θ E , which can be used to constrain the physical lens parameters. However, the microlens parallax vector, π E = (π rel /θ E )/(µ/µ), which is another observable constraining the physical lens parameters, cannot be securely measured because the event time scale, t E ∼ 19 days, is not long enough to produce detectable deviations induced by the orbital motion of Earth around the Sun (Gould 1992). Here µ represents the vector of the relative lens-source proper motion.

KMT-2021-BLG-0320
The lensing event KMT-2021-BLG-0320 was found on 2021 April 9 (HJD ′ ∼ 9313), when the event had not yet reached its peak, with the employment of the KMTNet AlertFinder system (Kim et al. 2018), which began full operation since the 2019 season. Two days after the detection, the event reached its peak with a magnification of A max ∼ 170, and then gradually declined until it reached its baseline of I base = 20.08. The source of the event lies in the two KMTNet prime fields of BLG01 and BLG41, toward which observations were conducted with a 30 min cadence for each field, and thus with a 15 min combined cadence. The areas covered by the two fields overlap except for about 15% of the area of each field filling the gaps between the chips of the camera. Because the data were taken from two fields of three telescopes, there are 6 data sets: KMTA (BLG01), KMTA (BLG41), KMTC (BLG01), KMT (BLG41), KMTS (BLG01), and KMTS (BLG41).
In Figure 4, we present the light curve constructed from the combination of the 6 KMTNet data sets. Although it would be difficult to notice the anomaly from a glimpse, we inspected the  -light curve because the event reached a very high magnification at the peak, near which the light curve is susceptible to perturbations induced by a planet (Griest & Safizadeh 1998). From this inspection, it was found that the light curve exhibited an anomaly that lasted for about 4 hours with a negative deviation with respect to a 1L1S model. See the enlarged view around the peak region of the light curve presented in the top panel of Figure 4. It is known that a planetary companion to a lens can induce anomalies with both positive and negative deviations, while a faint companion to a source can induce anomalies with only positive deviations (Gaudi 1998). We, therefore, modeled the light curve under the 2L1S interpretation. It is found that the anomaly is well described by a 2L1S model, in which the mass ratio between the binary lens components is very low. We found two sets of solutions with (s, q) close ∼ (0.77, 3.0×10 −4 ) and (s, q) wide ∼ (1.27, 2.9×10 −4 ). We designate the individual solutions as "close" and "wide" solutions, because the former solution has s < 1.0 and the latter has s > 1.0. In Table 4, we list the full lensing parameters for the two sets of solutions. The fits of the two solutions are nearly identical with χ 2 values that are equal to the first digit after the decimal point, indicating that the degeneracy between the two solutions is very severe. The fact that the binary separations of the two solutions follow the relation of s close × s wide ≃ 1.0 indicates that the sim-A&A proofs: manuscript no. ms ilarity between the solution stems from the close-wide degeneracy. It is known that the relation between the planet separations of the solutions under the offset degeneracy in Equation (1) applies to more general cases, including resonant case, than the s close × s wide ≃ 1.0 relation of the close-wide degeneracy, and we confirm this. In addition to the relation between s in and s out , Hwang et al. (2022) provided analytic formulas for the heuristic estimation of the source trajectory angle and the mass ratio; where τ anom denotes the duration of the planet-induced anomaly. We confirm that the heuristically estimated values of α and q match well those estimated from the modeling. It was found that the normalized source radius could not be constrained, not even the upper limit, and thus the line for ρ in Table 4 is left blank. Figure 5 shows the lensing configuration of KMT-2021-BLG-0320 for the close (upper panel) and wide (lower panel) solutions. It shows that the source passed the back-end region of the tiny central caustic induced by a planet, and this generated the negative deviation of the observed anomaly. The central caustics of the close and wide solutions are very similar, resulting in nearly identical deviations. Because the source passed well outside the caustic, finite-source effects could not be detected, and thus we do not mark an orange circle representing the source size. Microlens-parallax effects could not be securely detected because the event time scale, t E ∼ 13 days, is much shorter than the orbital period of Earth, that is, 1 year.

KMT-2021-BLG-1303 (MOA-2021-BLG-182)
The event KMT-2021-BLG-1303 was first found on 2021 June 14 (HJD ′ ∼ 9379) by the KMTNet survey, and later identified by the MOA survey on 2021 June 17 (HJD ′ ∼ 9382). The MOA survey designated the event as MOA-2021-BLG-182. A day after the first discovery, the event displayed an anomaly that lasted  for about 2 days. After the anomaly, the event reached peak at HJD ′ ∼ 9385.1 with a magnification of A max ∼ 46, and then gradually declined to the baseline of I base = 19.67. The light curve of KMT-2021-BLG-1303 constructed with the combined KMTNet and MOA data is shown in Figure 6. Compared to the previous two events, the light curve of this event displays a very obvious anomaly with a maximum deviation of ∆I ∼ 1.4 mag from a 1L1S model. The rapid brightening at HJD ′ ∼ 9381.0 indicates that the event experienced a caustic crossing. Because caustics form due to the multiplicity of lens components, we rule out the 1L2S origin of the anomaly and test only the 2L1S interpretation.
We find that the anomaly is well explained by a 2L1S model with a planetary companion lying close to the Einstein ring of the host. We find a unique solution without any degeneracy, and the binary parameters of the solution are (s, q) ∼ (1.029, 6.4 × 10 −4 ). Here again the planet-to-host mass ratio is of the order of 10 −4 like the two previous events. We double checked the uniqueness of the solution by thoroughly inspecting the region around the s and q parameters predicted by the relations in Equations (1) and (2), and confirmed that there is only a single solution. The full lensing parameters of the event are listed in Table 5.
The lensing configuration of the event is presented in Figure 7, which shows that the planet lies very close to the Ein- 6.98 ± 1.18 7.00 ± 1.28 98.8 ± 24.2 t 0,2 (HJD ′ ) --9394.740 ± 0.002 stein ring and induces a single six-sided resonant caustic near the primary of the lens. The source entered the caustic at HJD ′ ∼ 9381.0, and this produced the sharp rise of the light curve at the corresponding time. According to the model, the source exited the caustic at HJD ′ ∼ 9384.2, but the caustic-crossing feature at this epoch is not obvious in the lensing light curve due to the combination of the weak caustic and the poor coverage of the region. The pattern of the anomaly between the caustic entrance and exit deviates from a typical "U"-shape pattern because the source passed along the fold of the caustic. The light curve during the caustic entrance of the source was well resolved by the KMTA data, and thus the normalized source radius, ρ = (0.95 ± 0.08) × 10 −3 , is tightly constrained. However, it was difficult to constrain the microlens parallax because the time scale, t E ∼ 25 days, is not long enough and the photometric precision of the faint source with I ∼ 21 is not high enough to detect subtle deviations induced by microlens-parallax effects.

KMT-2021-BLG-1554
This short event with a time scale of t E ∼ 5 days and a faint baseline magnitude of I base = 22.39 was detected on 2021 Jul 1 (HJD ′ ∼ 9396), about two days after the event peaked at t 0 ∼ 9394.7, by the KMTNet survey. The source lies in the prime fields of BLG01 and BLG41, and thus there are 6 data sets from the three KMTNet telescopes. Among these data sets, the data set of the BLG41 field acquired by the KMTS telescope was not used in the analysis due to its poor photometric quality caused by bad seeing during the lensing magnification. Although the other data set from KMTS and those from KMTA were included in the analysis, the results are mainly derived from the KMTC data sets, which cover the peak region of the light curve.
The light curve constructed with the available KMTNet data sets is shown in Figure 8. In the peak region, it exhibits an anomaly that lasted for about 3 hours from a 1L1S model with a positive deviation of ∆I ∼ 0.5 mag. The anomaly appears both in the BLG01 (two data points) and BLG41 (3 points) data sets obtained from KMTC observations, confirming that the signal is real. Because a positive anomaly can be produced by a binary companion to either a lens or a source, we test both 2L1S and 1L2S models.
The lensing configurations of the event according to the close and wide solutions are shown in the upper and lower panels of Figure 9, respectively. The figure shows that the anomaly was produced by the crossing of the source over the sharp tip of the central caustic induced by a planet lying near the Einstein ring of the planet host. The gap between the caustic entrance and exit is much smaller than the source size, and thus the individual caustic-crossing features do not show up in the anomaly and instead the anomaly appears as a single bump. The incidence angle of the source on the binary axis is nearly 90 • , and thus the features on the rising and falling sides of the anomaly are symmetric.

Source stars and angular Einstein radii
Among the four analyzed planetary events, the anomalies of the three events (KMT-2017-BLG-0752, KMT-2021-BLG-1303, and KMT-2021-BLG-1554 were affected by finite-source effects, and thus the normalized source radii can be measured. In this section, we estimate the angular Einstein radii for these events. Estimating θ E from a measured ρ value requires one to specify the source type, from which the angular source radius θ * is deduced and the Einstein radius is determined by Although the normalized source radius and thus θ E cannot be measured for the event KMT-2021-BLG-0320, we specify the source type for the sake of completeness. The specification of the source type is done by estimating the color and brightness of the source. We measure the I and V-band magnitudes of the source from the regression of the data processed using the pyDIA photometry code (Albrow 2017) with the variation of the lensing magnification. We then place the source on the instrumental color-magnitude diagram (CMD) of stars lying around the source constructed using the same photometry code. Following the method of Yoo et al. (2004), the instrumental source color and magnitude, (V − I, I), are then calibrated using the centroid of red giant clump (RGC), for which its reddening and extinction-corrected (de-reddened) color and magnitude, (V − I, I) RGC,0 , are known, on the CMD as a reference. Figure 10 shows the positions of the source stars (marked by blue dots) with respect to those of the RGC centroids (red dots) on the CMDs of the individual events. For OGLE-2017-BLG-1691, the source color could not be securely measured not only because the V-band data sparsely covered the light curve during the lensing magnification but also because the photometry quality of these data is low, although the I-band magnitude was relatively well measured. In this case, we estimated the source color as the median color of stars with the same I-band magnitude offset from the RGC centroid in the CMD constructed from the images of Baade's window taken from the observations using the Hubble Space Telescope (Holtzman et al. 1998). From the offsets in color and magnitudes between the source and RGC centroid, ∆(V − I, I), together with the known de-reddened color and magnitude of the RGC centroid, (V −I, I) RGC,0 (Bensby et al. 2013;Nataf et al. 2013), we estimate the calibrated source color and magnitude as (V − I, I) 0 = (V − I, I) RGC,0 + ∆(V − I, I). In Table 7, we list the source colors and magnitudes of the individual events estimated from this procedure along with the spectral types of the source stars. It is found that the source of OGLE-2017-BLG-1691 is a turnoff star or a subgiant of a G spectral type, while the source stars of the other events are main-sequence stars with spectral types ranging from G to K.
For each event, the angular source radius and Einstein radius were estimated from the measured color and magnitude of the source. For this, the measured V − I color was converted into V − K color using the color-color relation of Bessell & Brett (1988), θ * value was interpolated from the (V − K)-θ * relation of Kervella et al. (2004), and then the angular Einstein radius was estimated using the relation in Equation (3). With the measured Einstein radius together with the event time scale, the value of the relative lens-source proper motion was assessed as In Table 8, we list the estimated values of θ * , θ E , and µ for the individual events. The θ E and µ values for the event KMT-2021-BLG-0320 are left blank because finite-source effects in the lensing light curve were not detected and subsequently the values of ρ, θ E , and µ could not be measured. We note that the angular Einstein radius of KMT-2021-BLG-1554, θ E ∼ 0.10 mas, is substantially smaller than those of the other events, and this together with its short time scale, t E ∼ 5 days, suggests that the mass of the lens would be very low. KMT-2021-BLG-1554 is the seventh shortest microlensing event with a bound planet. See Table 3 of Ryu et al. (2021).

Physical parameters
In this section, we estimate the physical parameters of the planetary systems. For the unique constraint of the physical lens parameters, one must simultaneously measure the extra observables of θ E and π E , from which the mass M and distance to the lens, D L , are determined as Here κ = 4G/(c 2 AU), π S = AU/D S , and D S denotes the distance to the source (Gould 1992(Gould , 2000. The values of the microlens parallax could not be measured for any of the events, and thus the physical parameters cannot be unambiguously determined from the relation in Equation (5). However, one can still constrain M and D L using the other measured observables of t E and θ E , which are related to the physical lens parameters as where π rel = AU(D −1 L − D −1 S ) denotes the relative parallax between the lens and source. The event time scales were well measured for all events, and the Einstein radii were assessed for three of the four events. In order to estimate the physical lens parameters with the constraints provided by the measured observables of the individual events, we conduct Bayesian analyses using a Galactic model.
In the Bayesian analysis, we started by producing many lensing events from a Monte Carlo simulation conducted with the use of a Galactic model, which defines the matter-density and kinetic distributions, and mass function of Galactic objects. We adopted the Galactic model constructed by Jung et al. (2021). In this model, the density distributions of disk and bulge objects are described by the Robin et al. (2003) and Han & Gould (2003) models, respectively. The kinematic distribution of disk objects is based on the modified version of the Han & Gould (1995) model, in which the original version is modified to reconcile the changed density model, that is, the Robin et al. (2003) model. The bulge kinematic distribution is modeled based on proper motions of stars in the Gaia catalog (Gaia Collaboration 2016. For the details of the density and kinematic distributions, see Jung et al. (2021). In the Galactic model, the mass function is constructed with the adoption of the initial mass function and the present-day mass function of Chabrier (2003) for the bulge and disk lens populations, respectively.
In the second step, we computed the time scales and Einstein radii of the simulated events using the relations in Equation (6), and then constructed posterior distributions of M and D L for the simulated events with the t E and θ E values that are consistent with the observables of the individual lensing events. For the three events OGLE-2017-BLG-1691, KMT-2021-BLG-1303, and KMT-2021-BLG-1554, we use the constraints of both t E and θ E , and for KMT-2021-BLG-0320, we apply the constraint of only t E because θ E is not measured for the event. Figures 11 and 12 show the Bayesian posteriors of M and D L of the individual events, respectively. In Table 9, we summarize the estimated values of the host and planet masses, M host and M planet , distance, and projected host-planet separation, a ⊥ = sD L θ E , for all of which medians are presented as representative A&A proofs: manuscript no. ms  tributions is not significant, although the mean lens mass for the former event is slightly bigger than that of the latter event due to the longer time scale and the uncertainty is smaller due to the additional constrain of θ E . The estimated masses indicate that all of the discovered planets have sub-Jovian masses. The planet masses KMT-2021-BLG-0320L, KMT-2021-BLG-1303L, and KMT-2021-BLG-1554L correspond to ∼ 0.10, ∼ 0.38, and ∼ 0.12 times of the mass of the Jupiter, and the planet mass of OGLE-2017-BLG-1691L corresponds to that of Uranus. To be noted among the planet hosts is that the estimated mass of the planetary system KMT-2021-BLG-1554L, M host ∼ 0.08 M ⊙ , corresponds to the boundary between a star and a brown dwarf (BD). Together with the previously detected planetary systems with very low-mass hosts 2 , the discovery of the system demonstrates the microlens-ing capability of detecting planetary systems with very faint or substellar hosts. Besides KMT-2021-BLG-1554L, the host stars of the other planetary systems are low-mass stars with masses lying in the range of ∼ [0.3-0.6] M ⊙ . Under the approximation that the snow line distance scales with the host mass as a sl ∼ 2.7(M/M ⊙ ) AU, all discovered planets lie beyond the snow lines of the hosts regardless of the solutions, demonstrating the high microlensing sensitivity to cold planets. The planetary systems lie in the distance range of ∼ [6.3-7.6] kpc, demonstrating the usefulness of the microlensing method in detecting remote planets.

Summary and conclusion
We presented the analyses of four planetary microlensing events OGLE-2017-BLG-1691, KMT-2021-BLG-0320, KMT-2021-BLG-1303, and KMT-2021-BLG-1554. The events share a common characteristic that the planetary signals appeared as anomalies with very short durations, ranging from a few hours to a couple of days, and they were clearly detected solely by the combined data of the high-cadence lensing surveys without additional data from followup observations.
From the detailed analyses of the events, it was found that the signals were generated by planets with low planet-to-host mass ratios: three of the planetary systems with mass ratios of the order of 10 −4 and the other with a mass ratio slightly greater than 10 −3 . In the histogram of microlensing planets presented in Figure 1, we mark the positions of the four planets discovered in this work. The estimated masses indicated that all discovered planets have sub-Jovian masses, in which the planet masses of KMT-2021-BLG-0320Lb, KMT-2021-BLG-1303Lb, and KMT-2021-BLG-1554Lb correspond to ∼ 0.10, ∼ 0.38, and ∼ 0.12 times of the mass of the Jupiter, and the mass of OGLE-2017-BLG-1691Lb corresponds to that of the Uranus. It was found that the host of the planetary system KMT-2021-BLG-1554L has a mass at around the boundary between a star and a brown dwarf. Besides this system, it was found that the host stars of the other planetary systems are low-mass stars with masses in the range of ∼ [0.3-0.6] M ⊙ . The discoveries of the planets well demonstrate the capability of the current high-cadence microlensing surveys in detecting low-mass planets.