Shear thickening in suspensions of particles with dynamic brush layers

Hojin Kim,∗a,b Michael van der Naald,a,c Finn A. Braaten,a,c Thomas A.
Witten,a,c Stuart J. Rowan,b,d and Heinrich M. Jaegera,c

a James Franck Institute, The University of Chicago, Chicago, Illinois 60637,
USA
b Pritzker School of Molecular Engineering, The University of Chicago,
Chicago, Illinois 60637, USA
c Department of Physics, The University of Chicago, Chicago, Illinois 60637,
USA
d Department of Chemistry, The University of Chicago, Chicago, Illinois 60637,
USA
∗<EMAIL_ADDRESS>
Control of frictional interactions among liquid-suspended particles has led to
tunable, strikingly non-Newtonian rheology via the formation of strong flow
constraints as particles come into close proximity under shear. Typically,
these frictional interactions have been in the form of physical contact,
controllable via particle shape and surface roughness. We investigate a
different route, where molecular bridging between nearby particle surfaces
generates a controllable ‘sticky’ friction. This is achieved with surface-
functionalized colloidal particles capable of forming dynamic covalent bonds
with telechelic polymers that comprise the suspending fluid. At low shear
stress this results in particles coated with a uniform polymer brush layer.
Beyond an onset stress $\sigma^{\ast}$ the telechelic polymers become capable
of bridging and generate shear thickening. Over the size range investigated,
we find that the dynamic brush layer leads to dependence of $\sigma^{\ast}$ on
particle diameter that closely follows a power law with exponent -1.76. In the
shear thickening regime, we observe an enhanced dilation in measurements of
the first normal stress difference $N_{1}$ and reduction in the extrapolated
volume fraction required for jamming, both consistent with an effective
particle friction that increases with decreasing particle diameter. These
results are discussed in light of predictions for suspensions of hard spheres
and of polymer-grafted particles.

##  1 Introduction

Shear thickening fluids are materials with a viscosity that rises with
increasing shear stress or shear rate [1]. Dense suspensions of solid
particles suspended in a Newtonian liquid have emerged as the prototypical
fluid exhibiting strong shear thickening and have been studied extensively
over the last two decades. Understanding the flow of suspensions is critically
important in both industrial applications and natural processes. Frictional
interaction between particles has been implicated as the dominant contribution
to dramatic increases in viscosity during shear thickening, whereby applied
stress leads to a growing population of particles interacting frictionally and
resisting the applied shear.[2, 3, 4, 5]

Advances in the understanding of shear thickening in dense suspensions have
been facilitated primarily by studying spherical particles as a model system.
This is the most straightforward for theoretical calculations and
computational modeling as well as being relatively simple to synthesize and
measure experimentally. Spherical particle suspensions also have the virtue of
making it relatively easy to isolate how changes to particle level
characteristics such as particle roughness, particle size, and particle
surface chemistry impact shear thickening. For example, by tuning particle
roughness [6, 7, 8, 9], varying contact friction [10], or incorporating
enhanced lubrication hydrodynamics in simulations [11, 12] it has been shown
how the shear thickening response can be controlled. Another way to alter the
shear thickening response is to tune the particle size, which controls the
stress $\sigma^{\ast}$ required to make frictional contact and thus to induce
shear thickening. For hard sphere suspensions this stress is found to scale as
$\sigma^{\ast}=\frac{F_{0}}{d^{2}}$, where the effective applied force
$\sigma^{\ast}d^{2}$ has to exceed the effective stabilizing force $F_{0}$ at
contact[13, 14]. Whenever local frictional contacts exist between contacting
particles, i.e., whenever $\sigma>\sigma^{\ast}$, relative movement is
hindered and this acts as a stress-activated constraint for the associated
degree of freedom [13, 15].

Figure 1: (a) Schematic of silica particles (grey) whose surfaces are
functionalized with thiol groups (red triangles) suspended in
p-nitrobenzalcyanoacetamide (N-BCAm)-endcapped poly(propylene glycol) (blue
polygons). The polymers form dynamic brushes at the particle surface through
the dynamic thia-Michael reaction, which convert to bridging bonds when
sheared. (b-d) Viscosity $\eta$ as a function of applied shear stress $\sigma$
for thiol-functionalized silica particles suspended in N-BCAm-endcapped
poly(propylene glycol) for different particle volume fraction $\phi$. Particle
diameter $d$ is 620nm (SP6, b), 1210nm (SP12, c), and 1930nm (SP19, d).
Vertical red dashed lines refer to the average onset stress for shear
thickening. Shaded area in left-top corner indicates inaccessible range due to
the rheometer minimum rate limit $\dot{\gamma}<0.001$ s-1.

Typically, such stress-induced constraints have their origin in physical
contact friction. However, friction can also be created by chemical bridging
between suitably functionalized surfaces.[16] For example, hydrogen bonding
between the particles can enhance the effective friction to the point of
inducing shear jamming [17]. However, hydrogen bonds are weak and have a very
short lifetime. Comparatively new and less studied are suspensions in which
the chemistry of particles and suspending polymeric liquid has been designed
to enable bridging via dynamic covalent bonds that can be significantly
stronger and have longer lifetimes than non-covalent interactions. These
polymers have groups at both ends that bond to the particle surfaces [18, 19].
The effective friction induced by such dynamic covalent bridging has been
termed “sticky friction” because the constraint is generated by solvent
molecules that chemically attach to contacting particle surfaces in a
reversible manner. These dynamic bonds establish quickly and release more
slowly, resulting in a longer lifetime. In addition to strong shear thickening
such suspensions also can exhibit anti-thixotopic responses owing to the slow
release [18, 19]. The focus of the present work is, however, only on the
suspension behavior with increasing applied stress and investigates how
dynamic bridging affects the onset of shear thickening, enables a cross-over
from lubrication- to friction-dominated rheological behavior, and changes the
propensity for jamming. Our experiments show how bridging interactions provide
a means of systematically tuning the effective friction in sheared suspensions
and, in particular, make it possible to change this effective friction simply
by changing particle size.

##  2 Experiments

###  2.1 Materials preparation

Five sets of silica spheres with hydroxyl surface groups were purchased from
Fiber Optic Center (New Bedford, MA) with diameters in the range 300 nm
$<d<2000$ nm. These particles were surface-functionalized with thiol groups as
reported previously (see Supplemental Material) [18]. The size of thiol-
functionalized particles was measured as $d=340\pm 4$ (SP3), $620\pm 10$
(SP6), $870\pm 10$ (SP9), $1210\pm 40$ (SP12), and $1930\pm 90$ nm (SP19)
using dynamic light scattering.

The suspending medium for these spheres was a telechelic polymeric liquid
comprised of p-nitrobenzalcyanoacetamide (N-BCAm)-endcapped poly(propylene
glycol) (number-averaged molecular weight, $M_{\mathrm{n}}\approx 5300$ g/mol)
[18]. A dynamic brush layer is formed on the surface of the particles via the
catalyst free, room temperature dynamic thia-Michael reaction between the
thiol and the BCAm group. In contrast to a covalently attached brush layer,
the dynamic brush continually forms and releases due to the kinetics of the
thia-Michael bond at room temperature (Fig. 1a). The equilibrium constant
$K_{\mathrm{eq}}$ of this thia-Michael reaction was measured as
$K_{\mathrm{eq}}\approx 8000$ M-1.[18] From a thermodynamic perspective, we
note that, for all suspensions reported here, the stoichiometric amount of the
BCAm Michael-acceptor group is at least 5-fold larger than the surface thiol
(1–2 thiol-per-nm2, see the nuclear magnetic resonance spectroscopy study in
the Supplemental Material). In this concentration limit and $K_{\mathrm{eq}}$,
the telechelic polymers form dynamic bonds with the surface thiol, resulting
in the brush layer.

###  2.2 Rheology measurement

Suspensions with different volume fractions $\phi$ were prepared by measuring
masses of particles and solvents and converting to volumes using their
densities. The rheology over all accessible shear stresses, typically
$\sigma=10^{-1}$ Pa to $\sigma=10^{3}$ Pa, was measured using a rheometer
(Anton Paar, MCR 302 or MCR 702) in cone and plate geometry (25 mm diameter,
2∘ angle). Owing to the high yield stress, the measurements were performed in
rate-controlled mode at low stress prior to the shear thickening regime and
otherwise using stress control. For suspensions with the smallest particles,
$d=340$ nm, the rheology was measured in rate-controlled mode throughout,
given the particularly high yield stress. The onset stress for shear
thickening was determined from the minimum viscosity.

##  3 Results

Fig. 1b shows the rheology of SP6 at different volume fractions $\phi$. For
$\phi=0.30$, the suspension exhibits a mild thinning and thickens slightly at
larger stresses. The magnitude of shear thickening increases with $\phi$, as
seen in the slope $d\log{\eta}/d\log{\sigma}$. At $\phi=0.45$ the slope
reaches 1, which corresponds to discontinuous shear thickening when viscosity
is plotted as a function of shear rate. For all $\phi$, the suspension starts
to shear thicken at $\sigma^{\ast}\approx 10$ Pa (red dashed line).

A second shear-thinning regime appears at high stress. While similar behavior
has been found in suspensions of anisotropic[20], irregular[21], or soft[22]
particles, this second thinning in the dynamic suspensions studied here is
most likely attributable to the breaking of the bridging interactions at high
shear stress. Probing larger particles, SP12 (Fig. 1c) and SP19 (Fig. 1d),
analogous behavior is observed. However, the shear thickening behavior for
given $\phi$ becomes weaker with increasing $d$. Similarly, the low stress
yielding behavior is weaker for larger $d$, even at larger $\phi$.

Figure 2: Onset stress for shear thickening $\sigma^{\ast}$ in suspensions
with dynamic covalent bonds (blue diamonds). Also shown are data from prior
studies [13, 23, 24, 14, 25, 26, 27, 28] using suspension with particles
having permanently attached covalent brush layers (black symbols). Solid lines
are least-squares fits.

Physically, the onset stress $\sigma^{\ast}$ signals the stress level at which
the particles have been sheared into sufficiently close distance, $h^{\ast}$,
to start forming frictional contact. In Fig. 2 we plot the onset stress for
shear thickening, $\sigma^{\ast}$, as a function of particle diameter $d$.
Over the size range measured, for our particles (blue symbols) this onset
closely follows a power law $\sigma^{\ast}\propto d^{\alpha}$ with best fit
exponent $\alpha=-1.76$. Also shown in that figure are results from prior work
[13, 23, 24, 14, 25, 26, 27, 28] with particles where the brush layer was
attached with permanent covalent bonds (black symbols), and for which a best
fit exponent $\alpha\approx-1.9$ has been reported [13]. We believe that one
of the reasons for our data to exhibit very little scatter is that the dynamic
brush layers is formed by the exact same N-BCAm-endcapped polymers for all
particle sizes.

Figure 3: (a) First normal stress difference $N_{1}$ as a function of shear
stress $\sigma$ for suspensions of 1930 nm particles (SP19) with volume
fractions $\phi=0.3$ to 0.52. The value measured at the lowest applied shear
stress, $N_{1,0}$, has been subtracted. (b) $N_{1}-N_{1,0}$ is plotted for
suspensions of particles with diameter $d=620$ nm (SP6), 1210 nm (SP12), and
1930 nm (SP19) at fixed $\phi=0.45$.

As $\sigma$ increases, particles are necessarily pushed closer together, thus
promoting bridging by the telechelic macromonomers via dynamic bonds at
$\sigma>\sigma^{\ast}$. This generates constraints to relative particle
movement, similar to the constraints arising from a network of frictional
contacts. In the following, we further explore the dynamic-bond-induced shear
thickening state through the first normal stress difference
$N_{1}\equiv\sigma_{\mathrm{xx}}-\sigma_{\mathrm{zz}}$. Here
$\sigma_{\mathrm{xx}}$ and $\sigma_{\mathrm{zz}}$ refer to the stress tensor
components in the shear flow direction (x) and its gradient direction (z).
Prior work has suggested that the sign of $N_{1}$ distinguishes states where
particles interact via lubrication ($N_{1}<0$) or form system-spanning
networks of frictional contacts that lead to dilation of the overall sample
volume ($N_{1}>0$) [29, 30, 3]. In Fig. 3a, we plot $N_{1}-N_{1,0}$ for the
largest particles, 1930nm, at different volume fractions $0.3\leq\phi\leq
0.52$, where the initial value $N_{1,0}$ at the lowest applied shear stress
has been subtracted. For $\phi\leq 0.5$, we find $N_{1}<0$ even beyond the
stress $\sigma^{\ast}\approx 1.3$ Pa that indicates the onset of shear
thickening. We interpret this as indicating that system-spanning networks of
frictional contacts can only form for $\sigma>\sigma^{\ast}$. This changes
once $\phi\geq 0.52$ where now the sign of $N_{1}$ reverses from negative to
positive as $\sigma\geq\sigma^{\ast}$.

The effect on $N_{1}$ of changing particle size is shown in Fig. 3b for fixed
$\phi=0.45$. Strikingly, the sign of $N_{1}$ in the shear thickening regime
reverses from negative to positive for particles with $d<1210$ nm (see
Supplemental Material for other particle packing fractions).

Figure 4: Viscosity $\eta_{\mathrm{t}}$ at the thickened state versus volume
fraction $\phi$ for SP6 (black), SP12 (blue), and SP19 (red). Lines are least-
squares fits to $\eta=A(1-\phi/\phi_{\mathrm{J}})^{-n}$. The jamming volume
fraction $\phi_{\mathrm{J}}$ are $\phi_{\mathrm{J,SP6}}=0.48$,
$\phi_{\mathrm{J,SP12}}=0.51$, and $\phi_{\mathrm{J,SP19}}=0.54$.

The size-dependent magnitude of the effective friction seen in the behavior of
$N_{1}$ is also observed in the approach to jamming as the particle volume
fraction $\phi$ is increased. This can be seen in plots of the high-stress
viscosity $\eta_{\mathrm{t}}$, measured at the upper end of the shear
thickened state, as a function of particle volume fraction $\phi$, which
diverges at a friction-dependent jamming volume fraction
$\phi^{\mu}_{\mathrm{J}}$ [31]. The $\phi^{\mu}_{\mathrm{J}}$ divergence
shifts to a lower value with increased interparticle friction and can
therefore be taken as another indicator of the effect of ‘sticky’ friction due
to dynamic bridging. In Fig. 4 we plot data for the three particle systems
measured and fit them to a power-law relation,
$\eta=A(1-\phi/\phi_{\mathrm{J}})^{-n}$ [32, 33, 31] to extract
$\phi^{\mu}_{\mathrm{J}}$. The fits, with $n$ = 1.7-1.8 (see Supplemental
Material), give $\phi^{\mu}_{\mathrm{J}}$ = 0.54, 0.51, and 0.48 for SP19,
SP12 and SP6, respectively.

##  4 Discussion

For particles coated with a polymer brush layer we can expect that shearing
them into contact produces a balance between the applied stress and the
counteracting stress from the stiffness of the brush. Balancing the force
$\sigma^{\ast}d^{2}$ with the restoring force derived from gradient of the
interparticle potential for brush-layer stabilized particles [14] gives

$\sigma^{\ast}d^{2}=-\left.\frac{\partial U(r)}{\partial
r}\right|_{r=d+h^{\ast}}.$ (1)

Different interparticle potentials will produce different scalings
$\sigma^{\ast}\propto d^{\alpha}$ of the onset stress as a function of
particle diameter $d$. For example, stabilization by steric interactions with
distance-independent forces $\frac{\partial U(r)}{\partial r}$ = const. leads
to $\sigma^{\ast}\propto d^{-2}$. The data shown Fig. 2 by the black symbols
therefore has been interpreted as the result of behavior close to the hard
sphere limit. Another much-studied limit is that for compressed semidilute
polymer brush layers[34, 14] with polymer brush thickness $L$ and surface-to-
surface distance $h$ (see Fig. 1a). In this case the stabilizing potential
gradient between polymer brush layers was found consistent with
$-\frac{\partial U(r)}{\partial r}\sim d^{0.25}$. This gradient gives a
scaling $\sigma^{\ast}\propto d^{-1.75}$[14]. Intriguingly, this closely
matches our data for dynamic brushes. Notably, however, the power law exponent
-0.25 in this scaling argument requires semidilute, good solvent conditions
for the brush, which are difficult to justify in our system, where we have a
very different situation: the particles are dispersed in a melt of N-BCAm
endcapped polymers and these polymers also form the dynamic brush layer.

The size dependence of the sign reversal in $N_{1}$ (Fig. 3b) suggests that
the strength of the effective friction produced by the bridging linkers
increases with decreasing $d$. To rationalize this, we consider how particle-
particle interactions due to polymer linkages depend on the total energy
density associated with bridging. This energy density ($\rho_{\mathrm{B}}$)
depends on the energy per bridging bond $\epsilon_{\mathrm{B}}$, the
interparticle contact area $A\sim Ld$ for $d\gg L$, and the number of
particles per unit volume, which is proportional to $d^{-3}$, and it scales as
$\rho_{\mathrm{B}}\sim\epsilon_{\mathrm{B}}Ld/d^{3}=\epsilon_{\mathrm{B}}Ld^{-2}$.
As a result, although telechelic polymers on larger particles with less curved
surfaces have a stronger attractive pair potential, particularly in the regime
of $d\gg R_{\mathrm{g}}$ [35], the overall bridging energy density scales
inversely with particle size, $\rho_{\mathrm{B}}\sim d^{-2}$. Hence, per unit
volume smaller particles are more strongly interacting by bridging and behave
more frictional, i.e., they experience stronger constraints to relative motion
under shear. This then decreases the volume fraction required for the
formation of system-spanning frictional networks and thus for the observation
of behavior with $N_{1}>0$. We also see from Fig. 3 that in suspensions where
$N_{1}<0$ initially, there can be a sign reversal at larger shear stresses. We
take this as an indication that telechelic linkers make it possible to
establish extended network structures even under large applied shear. However,
the decrease in viscosity seen in Fig. 1 at these large stress levels suggests
that such networks structures deform more easily under shear than hard spheres
would.

The fits to the data in Fig. 4 make it possible to estimate the associated
effective friction coefficient $\mu$, using the simulation data by Singh et
al.[36]. Taking $\phi_{\mathrm{J}}^{\infty}\approx 0.47$ as obtained for
infinite sliding friction $\mu_{\mathrm{s}}=\infty$ at fixed rolling friction
$\mu_{\mathrm{r}}=0.3$[15], the resulting values are $\mu$ = 3.5, 0.9, and 0.5
for SP6, SP12, and SP19, respectively. The finding that the extracted
effective friction coefficient can be larger than unity again highlights that
we are dealing with ‘sticky’ friction due to chemical bridging, as opposed to
physical contact friction for which $\mu\leq 1$ can typically be expected. We
note that it is possible to attribute some of the decrease in
$\phi^{\mu}_{\mathrm{J}}$ with particle diameter $d$ also to the brush layer
effectively increasing the volume fraction of particles in the suspension.
Accounting for the brush layer would then shift $\phi$ to
$\phi_{\mathrm{eff}}=\phi(1+2L/d)^{3}$. However, this shift is too small to
fully explain the observed particle size dependence in
$\phi^{\mu}_{\mathrm{J}}$ (see Supplemental Material).

##  5 Summary and Conclusions

We investigated shear thickening for a class of dense suspensions with
particles designed to interact through chemical bridging interactions. In
these systems silica particles are surface-functionalized with thiol groups
and suspended in a N-BCAm-endcapped poly(propylene glycol). This combination
of particle surface and solvent chemistry allows for a dynamic brush layer to
be established at the particle surface. At low shear this brush layer
stabilizes the particles. As larger stress is applied to the suspension,
particles come into sufficiently close proximity that dynamic bridging
interactions are formed. These bridges introduce an effective interparticle
friction and promote shear thickening. We find that the onset stress of shear
thickening for suspensions with dynamic brush layers scales with particle size
as $d^{-1.76}$. It is intriguing that this lines up closely with the
theoretical model for covalently-attached brush stabilized particles, which
predicts $d^{-1.75}$ but applies to a rather different regime. The degree to
which this model can be generalized is currently unknown. Our results show
that bridging interactions can give rise to dilatancy and form system-spanning
networks, just as in shear thickening suspensions where the particles interact
sterically through direct contact friction. However, owing to the increase in
the overall area of contact for fixed brush thickness, and in contrast to
steric interactions, the effective bridging-induced friction increases as
particle size decreases. Consequently, also the onset of jamming is shifted to
smaller $\phi^{\mu}_{\mathrm{J}}$ for smaller particles.

Our findings demonstrate that dynamic bridging interactions provide a direct,
highly tunable means for controlling the effective particle friction. For
fixed particle and solvent chemistry, this friction can be changed
straightforwardly by changing the particle size and it can become as large as
or even exceed physical friction produced by particle surface roughness. While
we reported on one specific dynamic covalent chemistry, the results introduce
a versatile way to design the non-Newtonian rheology of shear thickening
fluids, and we expect that this approach can be extended further by tuning the
molecular structure of the suspending polymer medium, including the molecular
weight (length scale of chemical ‘sticky’ friction), persistence length
(rolling friction), and thermodynamic reaction equilibrium (effective friction
lifetime).

## Author Contributions

HK, SJR, and HMJ conceived the research project; HK and FAB performed the
experiments; HK, MvdN, TAW, SJR, and HMJ analyzed the data; All authors
participated in the discussion and writing the manuscript.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

We thank Norman Wagner for fruitful discussion. The authors acknowledge
support from the University of Chicago Materials Research Science and
Engineering Center (MRSEC), which is supported by the National Science
Foundation under award DMR-2011854. This work was partially supported by
National Science Foundation under award DMR-2104694. Parts of this work were
carried out at the Soft Matter Characterization Facility and at the MRSEC
Characterization Facility at the University of Chicago (award number
DMR-2011854).

## Notes and references

  * [1] Jan Mewis and Norman J Wagner.  Colloidal suspension rheology.  Cambridge University Press, 2012.
  * [2] Ryohei Seto, Romain Mari, Jeffrey F Morris, and Morton M Denn.  Discontinuous shear thickening of frictional hard-sphere suspensions.  Physical Review Letters, 111(21):218301, 2013.
  * [3] John R Royer, Daniel L Blair, and Steven D Hudson.  Rheological signature of frictional interactions in shear thickening suspensions.  Physical Review Letters, 116(18):188301, 2016.
  * [4] Neil YC Lin, Ben M Guy, Michiel Hermes, Chris Ness, Jin Sun, Wilson CK Poon, and Itai Cohen.  Hydrodynamic and contact contributions to continuous shear thickening in colloidal suspensions.  Physical Review Letters, 115(22):228304, 2015.
  * [5] Chiao-Peng Hsu, Joydeb Mandal, Shivaprakash N Ramakrishna, Nicholas D Spencer, and Lucio Isa.  Exploring the roles of roughness, friction and adhesion in discontinuous shear thickening by means of thermo-responsive particles.  Nature Communications, 12(1):1477, 2021.
  * [6] Didier Lootens, Henri Van Damme, Yacine Hémar, and Pascal Hébraud.  Dilatant flow of concentrated suspensions of rough particles.  Physical Review Letters, 95(26):268302, 2005.
  * [7] Lilian C Hsiao, Safa Jamali, Emmanouil Glynos, Peter F Green, Ronald G Larson, and Michael J Solomon.  Rheological state diagrams for rough colloids in shear flow.  Physical Review Letters, 119(15):158001, 2017.
  * [8] Philippe Bourrianne, Vincent Niggel, Gatien Polly, Thibaut Divoux, and Gareth H McKinley.  Tuning the shear thickening of suspensions through surface roughness and physico-chemical interactions.  Physical Review Research, 4(3):033062, 2022.
  * [9] Chiao-Peng Hsu, Shivaprakash N Ramakrishna, Michele Zanini, Nicholas D Spencer, and Lucio Isa.  Roughness-dependent tribology effects on discontinuous shear thickening.  Proceedings of the National Academy of Sciences, 115(20):5117–5122, 2018.
  * [10] Yu-Fan Lee, Yimin Luo, Tianyi Bai, Carlos Velez, Scott C Brown, and Norman J Wagner.  Microstructure and rheology of shear-thickening colloidal suspensions with varying interparticle friction: Comparison of experiment with theory and simulation models.  Physics of Fluids, 33(3), 2021.
  * [11] Safa Jamali and John F Brady.  Alternative frictional model for discontinuous shear thickening of dense suspensions: Hydrodynamics.  Physical Review Letters, 123(13):138002, 2019.
  * [12] Mu Wang, Safa Jamali, and John F Brady.  A hydrodynamic model for discontinuous shear-thickening in dense suspensions.  Journal of Rheology, 64(2):379–394, 2020.
  * [13] Ben M Guy, Michiel Hermes, and Wilson CK Poon.  Towards a unified description of the rheology of hard-particle suspensions.  Physical Review Letters, 115(8):088304, 2015.
  * [14] Lakshmi-Narasimhan Krishnamurthy, Norman J Wagner, and Jan Mewis.  Shear thickening in polymer stabilized colloidal dispersions.  Journal of Rheology, 49(6):1347–1360, 2005.
  * [15] Abhinendra Singh, Grayson L Jackson, Michael van der Naald, Juan J de Pablo, and Heinrich M Jaeger.  Stress-activated constraints in dense suspension rheology.  Physical Review Fluids, 7(5):054302, 2022.
  * [16] Uri Raviv, Rafael Tadmor, and Jacob Klein.  Shear and frictional interactions between adsorbed polymer layers in a good solvent.  Journal of Physical Chemistry B, 105(34):8125–8134, 2001.
  * [17] Nicole M James, Endao Han, Ricardo Arturo Lopez de la Cruz, Justin Jureller, and Heinrich M Jaeger.  Interparticle hydrogen bonding can elicit shear jamming in dense suspensions.  Nature Materials, 17(11):965–970, 2018.
  * [18] Grayson L Jackson, Joseph M Dennis, Neil D Dolinski, Michael van der Naald, Hojin Kim, Christopher Eom, Stuart J Rowan, and Heinrich M Jaeger.  Designing stress-adaptive dense suspensions using dynamic covalent chemistry.  Macromolecules, 55(15):6453–6461, 2022.
  * [19] Hojin Kim, Mike van der Naald, Neil D. Dolinski, Stuart J. Rowan, and Heinrich M. Jaeger.  Dynamic-bond-induced sticky friction tailors non-newtonian rheology.  Soft Matter, 19:6797–6804, 2023.
  * [20] Eric Brown, Nicole A Forman, Carlos S Orellana, Hanjun Zhang, Benjamin W Maynor, Douglas E Betts, Joseph M DeSimone, and Heinrich M Jaeger.  Generality of shear thickening in dense suspensions.  Nature Materials, 9(3):220–224, 2010.
  * [21] Srinivasa R Raghavan and Saad A Khan.  Shear-thickening response of fumed silica suspensions under steady and oscillatory shear.  Journal of Colloid and Interface Science, 185(1):57–67, 1997.
  * [22] Safa Jamali, Arman Boromand, Norman Wagner, and Joao Maia.  Microstructure and rheology of soft to rigid shear-thickening colloidal suspensions.  Journal of Rheology, 59(6):1377–1395, 2015.
  * [23] Jan Mewis and Jan Vermant.  Rheology of sterically stabilized dispersions and latices.  Progress in Organic Coatings, 40(1-4):111–117, 2000.
  * [24] William J Frith, P d’Haene, R Buscall, and Joannes Mewis.  Shear thickening in model suspensions of sterically stabilized particles.  Journal of Rheology, 40(4):531–548, 1996.
  * [25] TA Strivens.  The shear thickening effect in concentrated dispersion systems.  Journal of Colloid and Interface Science, 57(3):476–487, 1976.
  * [26] MI Smith, R Besseling, ME Cates, and V Bertola.  Dilatancy in the flow and fracture of stretched colloidal suspensions.  Nature Communications, 1(1):114, 2010.
  * [27] L Isa.  Capillary flow of dense colloidal suspensions.  PhD thesis, University of Edinburgh, 2008.
  * [28] P d’Haene, Joannes Mewis, and GG Fuller.  Scattering dichroism measurements of flow-induced structure of a shear thickening suspension.  Journal of Colloid and Interface Science, 156(2):350–358, 1993.
  * [29] David R Foss and John F Brady.  Structure, diffusion and rheology of brownian suspensions by stokesian dynamics simulation.  Journal of Fluid Mechanics, 407:167–200, 2000.
  * [30] Nayoung Park, Vikram Rathee, Daniel L Blair, and Jacinta C Conrad.  Contact networks enhance shear thickening in attractive colloid-polymer mixtures.  Physical Review Letters, 122(22):228003, 2019.
  * [31] Romain Mari, Ryohei Seto, Jeffrey F Morris, and Morton M Denn.  Shear thickening, frictionless and frictional rheologies in non-brownian suspensions.  Journal of Rheology, 58(6):1693–1724, 2014.
  * [32] Bruno Andreotti, Jean-Louis Barrat, and Claus Heussinger.  Shear flow of non-brownian suspensions close to jamming.  Physical Review Letters, 109(10):105901, 2012.
  * [33] Edan Lerner, Gustavo Düring, and Matthieu Wyart.  A unified framework for non-brownian suspension flows and soft amorphous solids.  Proceedings of the National Academy of Sciences, 109(13):4798–4803, 2012.
  * [34] Glenn H Fredrickson and P Pincus.  Drainage of compressed polymer layers: dynamics of a” squeezed sponge”.  Langmuir, 7(4):786–795, 1991.
  * [35] Xiao-Xia Meng and William B Russel.  Telechelic associative polymers: Interaction potential and high frequency modulus.  Journal of Rheology, 50(2):169–187, 2006.
  * [36] Abhinendra Singh, Romain Mari, Morton M Denn, and Jeffrey F Morris.  A constitutive model for simple shear of dense frictional suspensions.  Journal of Rheology, 62(2):457–468, 2018.

# Electronic supplementary information:
Shear thickening in suspensions of particles with dynamic brush layers

Hojin Kim∗a,b, Michael van der Naalda,c, Finn A. Braatena,c,
Thomas A. Wittena,c, Stuart J. Rowanb,d, and Heinrich M. Jaegera,c

( aJames Franck Institute, The University of Chicago, Chicago, Illinois 60637,
USA
bPritzker School of Molecular Engineering, The University of Chicago, Chicago,
Illinois 60637, USA
cDepartment of Physics, The University of Chicago, Chicago, Illinois 60637,
USA
d Department of Chemistry, The University of Chicago, Chicago, Illinois 60637,
USA
∗<EMAIL_ADDRESS>
)

##  S1. Particle synthesis and characterization

Silica particles were purchased from Fiber Optic Center (New Bedford, MA) with
diameters 300 nm, 500 nm, 1 $\mathrm{\mu m}$, and 2 $\mathrm{\mu m}$. The
particle surface was functionalized by the procedure reported previously[1].
First, 10 g of silica particles were added to 400 mL of toluene. The
suspension was sonicated for 1 hr and stirred for another 30 min to homogenize
the suspension. Then 3-mercaptopropyl trimethoxysilane (MPTMS) was added to
the solution. The concentration of MPTMS was kept 10 MPTMS molecules per nm2
surface area of the added silica particles. The solution was then heated to
reflux and left for 24 hr. Toluene from the final solution was completely
removed using a rotary evaporator. The particle powder was washed with ethanol
three times by repeating sonication and centrifugation. The final particles
were left under vacuum for 24 hr.

Table S1:  Estimated thiol density of the particle surface. diameter |  340 nm |  620 nm |  870 nm |  1210 nm |  1930 nm
---|---|---|---|---|---
Thiol density of the particle surface |  1.09 SH/nm2 |  2.15 SH/nm2 |  2.19 SH/nm2 |  1.20 SH/nm2 |  1.31 SH/nm2

The surface density of thiol functional groups was estimated using nuclear
magnetic resonance (NMR) studies[1] using 1,3,5-trioxane as an internal
standard. First, 20 mg of the silica particles was added to 0.5 M NaOD/D2O
solution. The solution was stirred at 85 ∘C overnight to fully dissolve the
silica particle. The thiol density of the particle surface
$\rho_{\mathrm{SH}}$ was estimated by NMR. The peak intensity of the
(ONa)3SiCH2CH2CH2SH ($\delta=2.43$ ppm) was compared to that of the internal
standard (Table S1). The total number of the surface thiol $N_{\mathrm{SH}}$
at the particle with diameter $d$ in a suspension with packing fraction $\phi$
is estimated by

$N_{\mathrm{SH}}=\phi\frac{\pi d^{2}}{\pi d^{3}/6}\rho_{\mathrm{SH}},$

and the total number of the Michael-acceptor group is estimated by

$N_{\mathrm{MA}}=(1-\phi)\rho_{\mathrm{PPG}}N_{\mathrm{A}}n_{\mathrm{MA}}/M_{\mathrm{PPG}}$

with the Avogadro constant $N_{\mathrm{A}}$, the density
($\rho_{\mathrm{PPG}}\approx 1$ g/ml) and molecular weight
($M_{\mathrm{PPG}}=5300$ g/mol) of (N-BCAm)-endcapped poly(propylene glycol),
and the number of Michael-acceptor group $n_{\mathrm{MA}}=2$. The ratio of
Michael-acceptor to thiol $N_{\mathrm{MA}}/N_{\mathrm{SH}}$ in suspensions
with various particle volume fractions is plotted in Figure S1.

Figure S1: Stoichiometric ratio of Michael-acceptor to thiol in suspensions of different particle sizes is plotted over the particle volume fraction range studied in this study. Table S2: Parameters from a power-law relation, $\eta=A(1-\phi/\phi_{\mathrm{J}})^{-n}$. diameter |  $A$ |  $n$ |  $\phi_{\mathrm{J}}$
---|---|---|---
620 nm |  0.79 |  1.79 |  0.48
1210 nm |  0.96 |  1.69 |  0.51
1930 nm |  0.74 |  1.73 |  0.54
Figure S2: First normal stress difference $N_{1}$ for particles ($d=620$ nm,
SP6). Figure S3: First normal stress difference $N_{1}$ for particles
($d=1210$ nm, SP12). Figure S4: Time-dependent measurement of the first normal
stress difference $N_{1}$ for the suspension of particles with diameter
$d=620$ nm at $\phi=0.45$. Measurements were performed at each shear stress
after the consistent preshear and relaxation steps. Blue lines indicate the
$N_{1}$ at each shear stress from the stress-sweep measurement (30 s per each
stress). Figure S5: Viscosity $\eta_{\mathrm{t}}$ at the thickened state
versus the effective volume fraction $\phi_{\mathrm{eff}}=\phi(1+2L/d)^{3}$
for SP6 (black), SP12 (blue), and SP19 (red). The thickness of the brush layer
is set to half of the contour length $L_{\mathrm{c}}\approx 10$ nm. Lines are
least-squares fits to $\eta=A(1-\phi_{\mathrm{eff}}/\phi_{\mathrm{J}})^{-n}$.
The estimated jamming volume fraction $\phi_{\mathrm{J}}$ is
$\phi_{\mathrm{J,SP6}}=0.50$, $\phi_{\mathrm{J,SP12}}=0.52$, and
$\phi_{\mathrm{J,SP19}}=0.55$.

## Notes and references

  * [1] Carina IC Crucho, Carlos Baleizão, and José Paulo S Farinha.  Functional group coverage and conversion quantification in nanostructured silica by 1H NMR.  Analytical Chemistry, 89(1):681–687, 2017.


# KMT-2019-BLG-1715: planetary microlensing event with three lens masses and
two source stars

Cheongho Han1, Andrzej Udalski2, Doeon Kim1, Youn Kil Jung3, Chung-Uk Lee3,
Ian A. Bond4
(Leading authors),
and
Michael D. Albrow0007, Sun-Ju Chung3,7, Andrew Gould8,9, Kyu-Ha Hwang3, Hyoun-
Woo Kim3, Yoon-Hyun Ryu3, In-Gu Shin3, Yossi Shvartzvald10, Weicheng Zang11,
Jennifer C. Yee4, Sang-Mok Cha3,12, Dong-Jin Kim3, Seung-Lee Kim3,7, Dong-Joo
Lee3, Yongseok Lee3,12, Byeong-Gon Park3,7, Richard W. Pogge9, Chun-Hwey
Kim13, Woong-Tae Kim14
(The KMTNet Collaboration),
Przemek Mróz2,15, Michał K. Szymański2, Jan Skowron2, Radek Poleski2, Igor
Soszyński2, Paweł Pietrukowicz2, Szymon Kozłowski2, Krzysztof Ulaczyk16,
Krzysztof A. Rybicki2, Patryk Iwanek2, Marcin Wrona2, Mariusz Gromadzki2
(The OGLE Collaboration)
Fumio Abe17, Richard Barry18, David P. Bennett18,19, Aparna Bhattacharya18,19,
Martin Donachie20, Hirosane Fujii17, Akihiko Fukui21,22, Yoshitaka Itow17,
Yuki Hirao23, Rintaro Kirikawa23, Iona Kondo23, Man Cheung Alex Li20, Yutaka
Matsubara16, Yasushi Muraki16, Shota Miyazaki23, Clément Ranc18, Nicholas J.
Rattenbury20, Yuki Satoh23, Hikaru Shoji23, Haruno Suematsu23, Takahiro
Sumi23, Daisuke Suzuki23, Yuzuru Tanaka23, Paul J. Tristram24, Takeharu
Yamakawa16, Tsubasa Yamawaki23, Atsunori Yonehara25,
(The MOA Collaboration)
1 Department of Physics, Chungbuk National University, Cheongju 28644,
Republic of Korea,  2 Astronomical Observatory, University of Warsaw, Al.
Ujazdowskie 4, 00-478 Warszawa, Poland  3 Korea Astronomy and Space Science
Institute, Daejon 34055, Republic of Korea  4 Institute of Natural and
Mathematical Sciences, Massey University, Auckland 0745, New Zealand  5 Center
for Astrophysics $|$ Harvard & Smithsonian 60 Garden St., Cambridge, MA 02138,
USA  6 University of Canterbury, Department of Physics and Astronomy, Private
Bag 4800, Christchurch 8020, New Zealand  7 Korea University of Science and
Technology, 217 Gajeong-ro, Yuseong-gu, Daejeon, 34113, Republic of Korea  8
Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany
9 Department of Astronomy, The Ohio State University, 140 W. 18th Ave.,
Columbus, OH 43210, USA  10 Department of Particle Physics and Astrophysics,
Weizmann Institute of Science, Rehovot 76100, Israel  11 Department of
Astronomy and Tsinghua Centre for Astrophysics, Tsinghua University, Beijing
100084, China  12 School of Space Research, Kyung Hee University, Yongin,
Kyeonggi 17104, Republic of Korea  13 Department of Astronomy & Space Science,
Chungbuk National University, Cheongju 28644, Republic of Korea  14 Department
of Physics & Astronomy, Seoul National University, Seoul 08826, Republic of
Korea  15 Division of Physics, Mathematics, and Astronomy, California
Institute of Technology, Pasadena, CA 91125, USA  16 Department of Physics,
University of Warwick, Gibbet Hill Road, Coventry, CV4 7AL, UK  17 Institute
for Space-Earth Environmental Research, Nagoya University, Nagoya 464-8601,
Japan  18 Code 667, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
19 Department of Astronomy, University of Maryland, College Park, MD 20742,
USA  20 Department of Physics, University of Auckland, Private Bag 92019,
Auckland, New Zealand  21 Instituto de Astrofísica de Canarias, Vía Láctea
s/n, E-38205 La Laguna, Tenerife, Spain  22 Department of Earth and Planetary
Science, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo,
Bunkyo-ku, Tokyo 113-0033, Japan  23 Department of Earth and Space Science,
Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
24 University of Canterbury Mt. John Observatory, P.O. Box 56, Lake Tekapo
8770, New Zealand  25 Department of Physics, Faculty of Science, Kyoto Sangyo
University, 603-8555 Kyoto, Japan

###### Abstract

We investigate the gravitational microlensing event KMT-2019-BLG-171, of which
light curve shows two short-term anomalies from a caustic-crossing binary-
lensing light curve: one with a large deviation and the other with a small
deviation. We identify five pairs of solutions, in which the anomalies are
explained by adding an extra lens or source component in addition to the base
binary-lens model. We resolve the degeneracies by applying a method, in which
the measured flux ratio between the first and second source stars is compared
with the flux ratio deduced from the ratio of the source radii. Applying this
method leaves a single pair of viable solutions, in both of which the major
anomaly is generated by a planetary-mass third body of the lens, and the minor
anomaly is generated by a faint second source. A Bayesian analysis indicates
that the lens comprises three masses: a planet-mass object with $\sim
2.6~{}M_{\rm J}$ and binary stars of K and M dwarfs lying in the galactic
disk. We point out the possibility that the lens is the blend, and this can be
verified by conducting high-resolution followup imaging for the resolution of
the lens from the source.

###### Subject headings:

Gravitational microlensing (672); Gravitational microlensing exoplanet
detection (2147)

##  1\. Introduction

During the first-generation experiment (Udalski et al., 1994a; Alcock et al,
1997), microlensing observations were carried out with about one day cadence.
By employing wide-field cameras mounted on multiple telescopes, the cadence of
microlensing observations has been dramatically shortened, and now it reaches
15 minutes for the fields of the highest stellar concentration. With the
shortened cadence, the event detection rate has greatly increased from several
dozens/yr during the first-generation experiments to the current rate of more
than 3000 events/yr.

Light curves of lensing events often deviate from that of a single lens and a
single source (1L1S) event, which produces a smooth and symmetric light curve.
The most common cause of the deviation is the binarity of the lens: 2L1S
events (Mao & Paczyński, 1991). Interpreting the light curves of 2L1S events
was a challenging task at the time when such events were first detected, e.g.,
OGLE No. 7 (Udalski et al., 1994b), because the methodology for the analysis
of these anomalous events had not yet been developed. With the development of
various methodologies, e.g., ray-shooting technique (Bond et al., 2002; Dong
et al., 2009; Bennett, 2010) and contour integration algorithm (Gould &
Gaucherel, 1997; Bozza et al., 2018) developed for finite magnification
computations, followed by the theoretical understanding of the binary lensing
physics, 2L1S events are routinely detected and analyzed as they progress (Ryu
et al., 2010; Bozza et al., 2012).

With the increasing number of densely covered events, one occasionally
confronts events having light curves that exhibit extra deviations from those
of 2L1S events. One important cause of such extra deviations is the existence
of an additional lens or source component. There exist 12 confirmed events,
for which at least four objects, including the lens and source components, are
needed for the interpretations of the light curves. See Table 1 of Han et al.
(2021). Among these events, 9 events were produced by triple lens systems,
3L1S events. For two major reasons, modeling a 3L1S event is difficult. The
first is the complexity of the lensing behavior (Daněk & Heyrovský, 2015,
2019) and the resulting difficulties in analyzing deviations. Another
important obstacle in analyzing these events is the degeneracy problem, which
makes it difficult to find a correct solution among those resulting in similar
light curves despite dramatically different interpretations of the lens
system. Identifying various types of degeneracies and investigating their
causes are important not only to correctly interpret the lens system but also
to investigate similar degeneracies in following analyses.

We investigate the microlensing event KMT-2019-BLG-1715 and present the
analysis result. The event displays a light curve with a complex pattern, in
which there exist two short-lasting anomalies deviating from a caustic-
crossing 2L1S light curve. We model the observed light curve under various
interpretations and present the results.

The organization of the paper for the presentation of the analysis is as
follows. We address the acquisition and reduction of data in Sect. 2. We
depict the pattern of the anomalies that deviate from a 2L1S lensing light
curve in Sect. 3. We present various models describing the anomalies in Sect.
4. We inspect the origins of the degeneracies among the identified solutions
in Section 5. We characterize the source, and measure the Einstein radius in
Sect. 6. We probe a method that can resolve the degeneracies in the lensing
models in Sect. 7. In Sect. 8, we describe a Bayesian analysis conducted to
characterize the lens, and present the physical quantifies. We summarize the
analysis and make a conclusion in Sect. 9.

Table 1Data and error readjustment factors Data set | $k$ |  $\sigma_{\rm min}$ (mag) | $N_{\rm data}$
---|---|---|---
OGLE | 1.225 | 0.030 | 450
KMTA (BLG03) | 1.404 | 0.020 | 1307
KMTA (BLG43) | 0.812 | 0.050 | 1195
KMTA (BLG03) | 1.517 | 0.010 | 1109
KMTA (BLG43) | 1.217 | 0.020 | 1424
KMTC (BLG03) | 1.382 | 0.010 | 981
KMTS (BLG43) | 1.187 | 0.020 | 1314
Figure 1.—  Lensing light curve of KMT-2019-BLG-1715. The model obtained from
a 2L1S interpretation (close model) is drawn over the data points. The two top
panels show the zoom-in views around the times of the source star’s entrance
into and exit from the caustic according to the 2L1S model. The two insets in
the lower panel show the zoomed-in views around the two anomalies from the
2L1S model occurred at ${\rm HJD}^{\prime}\sim 8691.5$ (major anomaly) and
$\sim 8694.5$ (minor anomaly).

##  2\. Observations and Data

The event KMT-2019-BLG-1715 was found from the observations of a source star
lying at $({\rm RA},{\rm DEC})_{\rm J2000}=(18:01:29.21,-28:46:37.7)$, that
correspond to $(l,b)=(1^{\circ}\hskip-2.0pt.898,-2^{\circ}\hskip-2.0pt.914)$.
The apparent source magnitude at the baseline was $I_{\rm base}=18.65$.

The magnification of the source flux was discovered from the Korea
Microlensing Telescope Network (KMTNet: Kim et al., 2016) survey on 2019 July
20 (${\rm HJD}^{\prime}\equiv{\rm HJD}-2450000\sim 8684.7$), at which the
source became brighter by about $0.25$ magnitude from the baseline. The KMTNet
survey utilizes three telescopes that lie around the world in Australia
(KMTA), Chile (KMTC), and South Africa (KMTS). The aperture of each telescope
is 1.6 meter, and the field of view of the camera mounted on the telescope is
$4~{}{\rm deg}^{2}$. The source lies in the BLG03 and BLG43 fields, which
overlap with a small offset, and thus the data are composed of six sets, with
two sets from each telescope.

The event was independently found by two other lensing surveys of the Optical
Gravitational Lensing Experiment (OGLE: Udalski et al., 2015) and the
Microlensing Observations in Astrophysics (MOA: Bond et al., 2001) on 2019
July 28 (${\rm HJD}^{\prime}\sim 8692.7$) and August 1 (${\rm
HJD}^{\prime}\sim 8696.7$), respectively. The OGLE and MOA surveys designated
the event as OGLE-2019-BLG-1190 and MOA-2019-BLG-352, respectively. The OGLE
survey uses the 1.3 m Warsaw Telescope located at the Las Campanas Observatory
in Chile, and the MOA survey utilizes the 1.8 m MOA-II Telescope located at
Mt. John Observatory in New Zealand.

Figure 2.—  Distribution of $\Delta\chi^{2}$ on the $\log s$–$\log q$ plane.
Red, yellow, green, cyan, blue, and purple colors are used to designate points
with $<1n\sigma$, $<2n\sigma$, $<3n\sigma$, $<4n\sigma$, $<5n\sigma$, and
$<6n\sigma$, where $n=10$.

The data were processed utilizing the photometry pipelines of the individual
groups: Albrow (2017), Udalski (2003), and Bond et al. (2001) for the KMTNet,
OGLE, and MOA surveys, respectively. We note that these pipelines commonly
utilize the difference imaging algorithm that was developed to optimize
photometry in very crowded fields (Alard & Lupton, 1998). Although the MOA
data cover the caustic crossing at ${\rm HJD}^{\prime}\sim 8697.9$, which is
also densely covered by the KMTA data sets, they do not cover the anomalies of
our major interest to be discussed below. Furthermore, the photometric
uncertainties of the data are substantially larger than those of the other
data sets, and thus the MOA data are not used for analysis to minimize their
effects on modeling. The error bars of the data used in the analysis were
rescaled from those estimated from the automatized pipelines, $\sigma_{0}$, by
$\sigma=k(\sigma_{\rm min}^{2}+\sigma_{0}^{2})^{1/2}$, where $k$ denotes a
rescaling factor to normalize $\chi^{2}/{\rm dof}$, where dof represents the
degree of freedom, to unity, and $\sigma_{\rm min}$ is the scatter of data
(Yee et al., 2012). In Table 1, we list $(k,\sigma_{\rm min},N_{\rm data})$
for the six data sets. Here $N_{\rm data}$ denotes the number of data in each
data set.

Table 2Lensing parameters of 2L1S model Parameter | Close | Wide
---|---|---
$\chi^{2}$ | 9329.3 | 9313.7
$t_{0}$ (${\rm HJD}^{\prime}$) | $8697.039\pm 0.022$ | $8696.971\pm 0.016$
$u_{0}$ | $0.123\pm 0.002$ | $0.091\pm 0.001$
$t_{\rm E}$ (days) | $30.94\pm 0.31$ | $39.04\pm 0.35$
$s$ | $0.756\pm 0.006$ | $2.028\pm 0.002$
$q$ | $0.141\pm 0.002$ | $0.310\pm 0.006$
$\alpha$ (rad) | $4.945\pm 0.005$ | $4.913\pm 0.005$
$\rho$ ($10^{-3}$) | $0.65\pm 0.01$ | $0.57\pm 0.01$
Figure 3.—  Lensing configurations of the two 2L1S models. The blue dots
indicate the lens positions, the red figures represent the caustics, and the
arrowed line represents the source motion. The Einstein ring is marked by a
dotted circle.

##  3\. 2L1S interpretation

The observed light curve of KMT-2019-BLG-1715 is shown in Figure 1. At first
glance, the two pronounced caustic-crossing spikes, at ${\rm HJD}^{\prime}\sim
8694.7$ ($t_{2}$) and $\sim 8697.9$ ($t_{3}$), suggests that the event is a
usual caustic-crossing 2L1S event. Fitting the light curve with a 2L1S model
yields a pair of solutions with $(s,q)\sim(0.76,0.14)$ (“close” model) and
$(s,q)\sim(2.03,0.31)$ (“wide” model), where $s$ and $q$ denote the projected
separation (scaled to the Einstein radius $\theta_{\rm E}$) and the mass ratio
between the binary lens components, with masses $M_{1}$ and $M_{2}$,
respectively. We carry out the 2L1S modeling in two steps, in which the binary
parameters $(s,q)$ are searched for via a dense grid approach in the first
step, and locals appearing in the $\Delta\chi^{2}$ map are refined in the
second step. The ranges of the grid parameters are $-1.0\leq\log s\leq 1.0$
and $-3.0\leq\log q\leq 0.0$ with 60 divisions. Figure 2 shows the
$\Delta\chi^{2}$ distribution obtained from the grid search and the locations
of the two locals, i.e., close and wide solutions.

Figure 4.—  Light curves at the baseline during three different time spans:
280 days (top panel), 80 days (middle panel), and 10 days (bottom panel). The
colors of data points are same as those in Fig. 1.

We list the full lensing parameters of the two 2L1S solutions in Table 2. In
the table, the parameters $(t_{0},u_{0},t_{\rm E})$ denote the closest source
approach time to a lens reference position, the reference-source separation at
$t_{0}$, and the timescale of the event, respectively. As a lens reference, we
choose the center of mass for a close binary ($s<1.0$), and the effective lens
position, defined by Di Stefano & Mao (1996) and An & Han (2002), for a wide
binary ($s>1.0$). Because the light curve exhibits pronounced caustic-crossing
spikes, finite-source effects (hereafter “finite effects”) are considered in
modeling by including $\rho=\theta_{*}/\theta_{\rm E}$ (normalized “source
radius”) as an additional parameter. Here $\theta_{*}$ denote the angular
source radius. Hereafter, we denote $\theta_{*}$ and $\theta_{\rm E}$ without
the notation of “angular”. In the computation of finite magnifications, we
take limb-darkening effects into consideration assuming that the brightness
profile on the source surface varies as $\Sigma_{\lambda}\propto
1-(1-3\cos\phi/2)\Gamma_{\lambda}$. Here $\Gamma_{\lambda}$ denotes the
coefficient of limb darkening, and $\phi$ denotes the angle between the two
lines extending from the source center, one toward an observer and the other
toward the source surface. Based on the spectral type, to be discussed in
Section 6, we adopt the $V$\- and $I$-band coefficients of
$(\Gamma_{V},\Gamma_{I})=(0.62,0.45)$. The two solutions (close and wide) are
degenerate, and the degeneracy is originated from the close–wide degeneracy
known by Dominik (1999). The lensing configurations, showing the trajectory of
the source relative to the lens components and caustic, are provided in Figure
3. The model curve of the 2L1S solution (for the close model) is shown in
Figure 1.

Figure 5.—  Zoomed-in views of the light curve around the major (at around
$t_{1}$) and minor (at around $t_{2}$) anomalies. For each epoch, the lower
five panels show the residuals from the identified five sets of degenerate
solutions: ELES (close), ESESin (wide), ESELout (wide), ESELin (wide), and
ESELout (wide) solutions. For each set, there are two solutions with
$s_{2}<1.0$ (close) and $s_{2}>1.0$ (wide), and the presented residuals are
for the solution yielding a better fit among the close–wide solutions. Drawn
over the data points are models of five degenerate solutions. The degeneracies
among the solutions are severe, resulting in similar model curves. The only
significant difference is the extra caustic entrance at $t_{2}$ corresponding
to the ESESin model.

Although the 2L1S model approximately describes the overall features of the
observed light curve, a close look reveals that the data exhibit two short-
lasting anomalies from the base 2L1S model. The first anomaly appears at ${\rm
HJD}^{\prime}\sim 8691.5$ ($t_{1}$) with $\Delta I\sim 0.7$ mag deviation from
the 2L1S model. The second anomaly appears just before the second caustic
crossing at $\sim 8694.5$ ($t_{2}$) with $\Delta I\sim 0.1$ mag deviation. We
refer to the former and latter anomalies as the “major anomaly” and “minor
anomaly”, respectively. The two insets in the lower panel of Figure 1 show the
zoom-in views of the anomalies. The major anomaly exhibits both rising and
falling parts during its very short duration of $\sim 3$ hours. The duration
of the minor anomaly is short as well, but precisely estimating the duration
is difficult because the coverage of the anomaly is incomplete.

We check the possibility that the extra anomalies arise due to photometric
artifacts or stellar variability. First, the possibility of the photometric
artifact, such as the change in transparency or the passage of a Solar System
object across the source, is very unlikely because both anomalies are covered
by multiple remotely separated telescopes and the data during the anomalies
show consistent patterns of variations. Second, the possibility of the source
variability is also unlikely because the stellar type of the source is a main
sequence, for which the chance of light variation is very small. We check the
source variability by inspecting the light curve at the baseline. Figure 4
shows the baseline light curve during three different time spans: 280 days
(top panel), 80 days (middle panel), and 10 days (bottom panel). It is found
that the light curve does not show any variation in all inspected time ranges.

Table 3Lensing parameters of ESES models Parameter | ESESin | ESESout
---|---|---
| Close | Wide | Close | Wide
$\chi^{2}$ | 7737.0 | 7696.1 | 7771.2 | 7701.1
$t_{0,1}$ (${\rm HJD}^{\prime}$) | $8697.047\pm 0.019$ | $8697.153\pm 0.017$ | $8697.046\pm 0.015$ | $8697.113\pm 0.016$
$u_{0,1}$ | $0.103\pm 0.002$ | $0.107\pm 0.001$ | $0.102\pm 0.001$ | $0.096\pm 0.001$
$t_{0,2}$ (${\rm HJD}^{\prime}$) | $8693.518\pm 0.035$ | $8694.716\pm 0.059$ | $8693.525\pm 0.029$ | $8694.400\pm 0.052$
$u_{0,2}$ | $0.205\pm 0.002$ | $0.368\pm 0.002$ | $0.204\pm 0.002$ | $0.304\pm 0.002$
$t_{0,3}$ (${\rm HJD}^{\prime}$) | $8696.383\pm 0.037$ | $8697.270\pm 0.058$ | $8696.669\pm 0.047$ | $8697.865\pm 0.105$
$u_{0,3}$ | $0.188\pm 0.003$ | $0.309\pm 0.003$ | $0.222\pm 0.003$ | $0.357\pm 0.007$
$t_{\rm E}$ (days) | $32.58\pm 0.28$ | $35.43\pm 0.17$ | $32.66\pm 0.05$ | $37.31\pm 0.25$
$s$ | $0.696\pm 0.004$ | $1.988\pm 0.001$ | $0.693\pm 0.003$ | $2.078\pm 0.003$
$q$ | $0.177\pm 0.003$ | $0.349\pm 0.004$ | $0.178\pm 0.002$ | $0.385\pm 0.007$
$\alpha$ (rad) | $4.997\pm 0.005$ | $4.950\pm 0.004$ | $4.999\pm 0.004$ | $4.957\pm 0.005$
$\rho_{1}$ ($10^{-3}$) | $0.59\pm 0.01$ | $0.63\pm 0.01$ | $0.60\pm 0.01$ | $0.58\pm 0.01$
$\rho_{2}$ ($10^{-3}$) | $0.47\pm 0.04$ | $0.45\pm 0.04$ | $0.42\pm 0.04$ | $0.43\pm 0.04$
$\rho_{3}$ ($10^{-3}$) | – | – | – | –
$q_{F,2}$ | $0.051\pm 0.002$ | $0.060\pm 0.002$ | $0.050\pm 0.001$ | $0.059\pm 0.002$
$q_{F,3}$ | $0.020\pm 0.003$ | $0.021\pm 0.003$ | $0.029\pm 0.004$ | $0.040\pm 0.006$

##  4\. Interpreting the anomalies

We investigate lensing models that can explain the observed anomalies from the
2L1S models. From this investigation, we find that adding a fourth body
(either a source or a lens) can explain one of the two anomalies, but not both
simultaneously. Solutions explaining both anomalies can only be found from
combinations of an extra lens (EL) and an extra source (ES), and we find five
sets of such solutions resulting from various types of degeneracy. For each
set, there exist two solutions with $s_{2}>1.0$ and $s_{2}<1.0$ resulting from
the close–wide binary degeneracy, and thus there are ten solutions in total.
The details of the individual models are discussed in the following
subsections.

For the designation of the individual solutions to be discussed, we use the
notations “ELES” (extra lens and extra source), “ESES” (extra source and extra
source), and “ESEL” (extra source and extra lens). In these notations, the
first “L” (or “S”) indicates that an extra lens (an extra source) is included
in the model to explain the major anomaly, while the second “L” (or “S”)
denotes that an extra lens (or source) is included to describe the minor
anomaly. For example, “ELES” denotes a model, in which the major and minor
anomalies are explained by adding an extra lens and an extra source to the
base 2L1S solution, respectively, and thus the model consists of three lens
components and two source stars.

Figure 6.—  Lensing configurations for the four ESES models: ESESin (close),
ESESin (wide), ESESout (close), and ESESout (wide). For each solution, the
panel on the left shows the wide view including the lens positions, and the
panel on the right shows the central region. The positions of the lens
components are marked by two filled dots ($M_{1}$ and $M_{2}$). The three
lines with arrows denote the motion of the primary ($S_{1}$, black), second
($S_{1}$, blue), and third ($S_{1}$, magenta) source stars.

###  4.1. ESES solutions

The ESES solutions explain both the major and minor anomalies with the
addition of two extra source stars to the base 2L1S model. In the modeling
considering additional source stars, the initial lensing parameters related to
the first source are adopted from those of the 2L1S model, and the initial
parameters related to the other source stars are assigned considering the
locations and magnitudes of the anomalies. According to these solutions, then,
there are two lens masses ($M_{1}$ and $M_{2}$) and three source stars
($S_{1}$, $S_{2}$, and $S_{3}$): 2L3S model. Besides the close-wide
degeneracy, the ESES interpretation additionally suffers from a degeneracy
arising due to the ambiguity in the trajectory of $S_{3}$, resulting in four
degenerate solutions in total. The two sets of the degenerate solutions
resulting from the latter type degeneracy, ESESin and ESESout, commonly
explain the origin of the major anomaly as the second source ($S_{2}$) star’s
caustic crossing over the tip of the binary-induced caustic. According to the
ESESin and ESESout solutions, the tertiary source ($S_{3}$) star, that is
introduced to explain the minor anomaly, passes the inner and outer region of
the caustic, respectively. The degeneracy between the two solutions results
from the incomplete coverage of the minor anomaly.

Table 4Lensing parameters of ESEL models Parameter | ESELin | ESELout
---|---|---
| Close | Wide | Close | Wide
$\chi^{2}$ | 7713.6 | 7668.5 | 7741.3 | 7667.9
$t_{0,1}$ (${\rm HJD}^{\prime}$) | $8697.094\pm 0.015$ | $8697.071\pm 0.016$ | $8697.141\pm 0.014$ | $8697.023\pm 0.015$
$u_{0,1}$ | $0.109\pm 0.002$ | $0.092\pm 0.001$ | $0.121\pm 0.001$ | $0.078\pm 0.002$
$t_{0,2}$ (${\rm HJD}^{\prime}$) | $8693.547\pm 0.034$ | $8694.239\pm 0.057$ | $8693.599\pm 0.033$ | $8693.982\pm 0.066$
$u_{0,2}$ | $0.216\pm 0.003$ | $0.292\pm 0.004$ | $0.232\pm 0.002$ | $0.229\pm 0.010$
$t_{\rm E}$ (days) | $32.05\pm 0.39$ | $38.35\pm 0.34$ | $30.57\pm 0.28$ | $42.92\pm 1.03$
$s_{2}$ | $0.716\pm 0.004$ | $2.089\pm 0.009$ | $0.739\pm 0.003$ | $2.227\pm 0.026$
$q_{2}$ | $0.164\pm 0.002$ | $0.372\pm 0.007$ | $0.162\pm 0.002$ | $0.398\pm 0.012$
$\alpha$ (rad) | $4.991\pm 0.005$ | $4.947\pm 0.005$ | $4.992\pm 0.005$ | $4.954\pm 0.004$
$s_{3}$ | $1.178\pm 0.006$ | $0.952\pm 0.004$ | $1.098\pm 0.003$ | $0.891\pm 0.005$
$q_{3}$ ($10^{-3}$) | $0.26\pm 0.05$ | $0.06\pm 0.01$ | $0.05\pm 0.01$ | $0.03\pm 0.01$
$\psi$ (rad) | $0.734\pm 0.004$ | $0.902\pm 0.003$ | $0.711\pm 0.003$ | $0.929\pm 0.005$
$\rho_{1}$ ($10^{-3}$) | $0.62\pm 0.01$ | $0.57\pm 0.01$ | $0.65\pm 0.01$ | $0.50\pm 0.02$
$\rho_{2}$ ($10^{-3}$) | $0.52\pm 0.04$ | $0.43\pm 0.04$ | $0.48\pm 0.04$ | $0.34\pm 0.04$
$q_{F}$ | $0.052\pm 0.002$ | $0.059\pm 0.002$ | $0.049\pm 0.002$ | $0.057\pm 0.002$

The lensing configurations of the four ESES solutions are presented in Figure
6, in which the trajectories of the second and third source stars are marked
by $S_{2}$, and $S_{3}$, respectively. The model curves and residuals from the
models in the regions around the anomalies are shown in Figure 5. We note that
the presented models are the wide solutions, which yield better fits than the
corresponding close solutions. The minor anomaly according to the ESESin
solution is characterized by a U-shape trough pattern because $S_{3}$ passes
the inner region of the caustic. The lensing parameters of the ESESin and
ESESout, including both the close and wide solutions, are listed in Table 3.
The flux ratios of the secondary and tertiary source stars to the primary
source are $q_{F,2}=F_{S,2}/F_{S,1}\sim 0.05$–0.06 and
$q_{F,3}=F_{S,3}/F_{S,1}\sim 0.02$–0.05, respectively. Here $F_{S,i}$ and
$(t_{0,i},u_{0,i},\rho_{i})$ indicate the flux and the lensing parameters
related to the $i$th source star, respectively. It is found that $\rho_{3}$
cannot be firmly determined because the coverage of the minor anomaly is
incomplete.

Figure 7.—  Lensing configurations of the four ESEL solutions: ESELin (close),
ESELin (wide), ESELout (close), and ESELout (wide). Notations are same as
those in Fig. 6, except that there are three lens components ($M_{1}$,
$M_{2}$, and $M_{3}$) and two source stars ($S_{1}$ and $S_{2}$).

###  4.2. ESEL solutions

According to the ESEL solutions, the lens and source of the event comprise
three lens masses ($M_{1}$, $M_{2}$, and $M_{3}$) and two source stars
($S_{1}$ and $S_{2}$), respectively: 3L2S model. Finding solutions is carried
out in two steps. In the first step, we fit the minor anomaly by adding an
extra lens component: 3L1S model. Because the overall pattern of the light
curve is described by a 2L1S model, the anomaly can be treated as a
perturbation to the 2L1S model (Bozza, 1999; Han et al., 2001). Under this
assumption, we search for a 3L1S model by fixing the parameters of the 2L1S
model ($s_{2}$, $q_{2}$, and $\alpha$), and then conducting a grid search for
the parameters describing the third lens component ($s_{3}$, $q_{3}$, and
$\psi$). Here ($s_{3}$, $q_{3}$) are the separation and mass ratio of $M_{3}$
with respect to $M_{1}$, respectively, and $\psi$ indicates the $M_{3}$
orientation angle, which is measured from the $M_{1}$–$M_{2}$ axis in a
clockwise sense centered at the position of $M_{1}$. The parameter ranges are
$-1.0\leq\log s_{3}\leq 1.0$, $-5.0\leq\log q_{3}\leq 1.0$, and
$0\leq\psi<2\pi$, and they are divided into 50, 50, and 180 grids,
respectively. The solutions are refined by gradually narrowing down the ranges
of the grid parameters, and then by releasing all parameters as free
parameters. In the second step, we fit the major anomaly by introducing an
extra source $S_{2}$ based on the 3L1S solution found in the first step.
Adding an additional source component requires us to include additional
parameters, including $(t_{0,2},u_{0,2},\rho_{2},q_{F})$, where $q_{F}$
represents the $S_{1}/S_{2}$ flux ratio.

Under the ESEL interpretation, we find four solutions, in which there exist a
pair of solutions for each of the close and wide solutions. For all of these
solutions, the major anomaly is generated by the crossing of $S_{2}$ over the
$M_{2}$-induced caustic, and the minor anomaly is generated by the passage of
$S_{1}$ through the region around the tiny $M_{3}$-induced caustic. Under this
interpretation, there are two solutions, designated as ESELin and ESELout
solutions, for each of the close and wide solutions. The difference between
the ESELin and ESELout solutions is that the first source ($S_{1}$) passes the
inner side (with respect to $M_{1}$) of the caustic induced by $M_{3}$ for the
ESELin solution, while $S_{1}$ passes the outer side of the caustic for the
ESELout solution. Figure 7 shows the lensing configurations of the four ESEL
solutions.

The model curves of the wide ESELin and ESELout solutions around the anomalies
are shown in Figure 5, and the parameters of the four ESEL solutions are
listed in Table 4. We note that the mass ratio between $M_{3}$ and $M_{1}$ of
the lens, $q_{3}\sim(4$–$26)\times 10^{-5}$, is very small, indicating that
$M_{3}$ is a very low-mass planet according to these solutions. Due to the
small mass ratio, the caustic induced by $M_{3}$ is much smaller than the
caustic induced by $M_{2}$.

Table 5Lensing parameters of ELES models Parameter | Close | Wide
---|---|---
$\chi^{2}$ | 7689.8 | 7713.7
$t_{0,1}$ (${\rm HJD}^{\prime}$) | $8696.722\pm 0.012$ | $8696.819\pm 0.009$
$u_{0,1}$ | $0.056\pm 0.001$ | $0.053\pm 0.001$
$t_{0,2}$ (${\rm HJD}^{\prime}$) | $8695.773\pm 0.027$ | $8696.341\pm 0.054$
$u_{0,2}$ | $0.147\pm 0.001$ | $0.181\pm 0.002$
$t_{\rm E}$ (days) | $43.96\pm 0.06$ | $51.93\pm 0.08$
$s_{2}$ | $0.551\pm 0.003$ | $2.524\pm 0.001$
$q_{2}$ | $0.246\pm 0.006$ | $0.461\pm 0.002$
$\alpha$ (rad) | $4.894\pm 0.005$ | $4.897\pm 0.005$
$s_{3}$ | $1.054\pm 0.001$ | $0.904\pm 0.001$
$q_{3}$ ($10^{-3}$) | $4.01\pm 0.16$ | $2.30\pm 0.10$
$\psi$ (rad) | $1.348\pm 0.008$ | $1.341\pm 0.006$
$\rho_{1}$ ($10^{-3}$) | $0.44\pm 0.01$ | $0.41\pm 0.01$
$\rho_{2}$ ($10^{-3}$) | – | –
$q_{F}$ | $0.043\pm 0.005$ | $0.042\pm 0.005$
Figure 8.—  $\Delta\chi^{2}$ map on the $s_{3}$–$q_{3}$ plane for the close
and wide 3L1S models. The colors of points follows the same coding as in Fig.
2.

###  4.3. ELES solutions

We find another 3L2S solutions, in which the way of interpreting the observed
anomalies is different from that of the 3L2S solutions described in the
previous subsection. We find these solutions by first fitting the major
anomaly from a 3L1S modeling, and then fitting the minor anomaly with the
introduction of an extras source $S_{2}$. We note that the ESEL solutions,
discussed in the previous subsection, are obtained by fitting the minor
anomaly first from a 3L1S modeling. According to the new 3L2S solutions, the
major and minor anomalies are explained with the addition of $M_{3}$ and
$S_{2}$ to the baseline 2L1S model, i.e., ELES models. For each of the close
and wide solutions, there is no additional degeneracy, and thus there are two
ELES solutions. In Figure 8, we present the $\Delta\chi^{2}$ map on the
$s_{3}$–$q_{3}$ plane.

We list the lensing parameters of the two ELES solutions (close and wide) in
Table 5. The estimated mass ratio, $q_{3}\equiv M_{3}/M_{1}\sim(2$–$4)\times
10^{-3}$, is very low, suggesting that $M_{3}$ is a planetary-mass object
belonging to the $M_{1}$–$M_{2}$ binary system according to the ELES model.
The separation and the orientation angle of the planet from $M_{1}$ is
$s_{3}\sim 1$ and $\psi=1.34$ (77.3∘), respectively. The value of $\rho_{2}$
is not presented in the table because $S_{2}$ does not involve with caustic
crossings, resulting in no finite effect.

Figure 9.—  Lensing configurations of the two ELES solutions: ELES (close) and
ELES (wide).

The model curve and the residuals from the ELES solutions at around the times
of the anomalies are shown in Figure 5. The lensing configurations are
presented in Figure 9. It is found that the configuration is greatly different
from that of ESEL solutions, although the lens systems of both ESEL and ELES
solutions constitute the same numbers of the lens and source components, i.e.,
3L2S model. We note that the caustics exhibit a self-intersecting pattern,
which is a characteristic pattern of a 3L system (Gaudi et al., 1998). From
the comparison of the lensing configurations with that of the 2L1S model
(shown in Figure 3), one finds that $M_{3}$ induces an additional set of
caustics, that overlaps with the one induced by $M_{2}$. According to the
solutions, the major anomaly arises due to the crossing of $S_{1}$ over one
tip of the caustic induced by $M_{3}$, while the minor anomaly is produced by
the passage of $S_{2}$ just outside of the tip of caustic induced by $M_{2}$.
The $I$-band flux ratio between $S_{1}$ and $S_{2}$ is $q_{F}\equiv
F_{S,2}/F_{S,1}\sim 0.042$–0.043, indicating that the second source is a faint
star.

##  5\. Origins of degeneracies

Despite the great differences in the interpretations, it is found that all the
degenerate solutions provide reasonably good fits to the data. In Table 6, we
compared the $\chi^{2}$ values of the ten models. It shows that the close
model yields a better fit than the wide model for the ELES solutions, while
the fits of wide models are better than the fits of the corresponding close
models for the other solutions. The middle panel of Figure 10 displays the
cumulative $\Delta\chi^{2}$ distributions from the base 2L1S model,
$\Delta\chi^{2}=\chi^{2}-\chi^{2}_{\rm 2L1S}$, for the five solutions, in
which we show the distribution of the model providing a better fit among each
pair of the close and wide solutions. The distributions show that all models
well describe both the major and minor anomalies, for which the times of the
anomalies, that is, $t_{1}$ and $t_{2}$, are marked by dotted lines. The
presented distributions of the solutions are very alike, making it difficult
to distinguish the distributions. For the better presentation of the
differences in the fits among the models, we additionally show the
distributions of the $\chi^{2}$ differences from the best-fit solution, i.e.,
$\Delta\chi^{2}=\chi^{2}-\chi^{2}_{\rm ESEL_{out}}$, in the bottom panel.

Table 6$\chi^{2}$ values of models

Solution |  $\chi^{2}$ ($\Delta\chi^{2}$)
---|---
| Close | Wide
ESESin | 7737.0 (69.1) | 7696.1 (28.2)
ESESout | 7771.2 (103.3) | 7701.1 (33.2)
ESELin | 7713.6 (45.7) | 7668.5 (0.6)
ESELout | 7741.3 (73.4) | 7667.9
ELES | 7689.8 (21.9) | 7713.7 (45.8)

Note. — The numbers in the parenthesis represent the $\chi^{2}$ difference
from the best-fit solution, i.e., $\Delta\chi^{2}=\chi^{2}-\chi^{2}_{\rm
ESEL_{\rm out}}$, where the best fit is given by the wide ESELout model.

With the multiple interpretations of the event, we probe the origins of the
degeneracy. As mentioned, the first type arises due to the ambiguity in the
separation between $M_{1}$ and $M_{2}$, i.e., the degeneracy between the close
and wide solutions (Griest & Safizadeh, 1998; Dominik, 1999; Albrow et al.,
2001). The second type degeneracy arises because the individual short-term
anomalies can be described either by an extra low-mass lens component or an
extra faint source component, i.e, the degeneracy between EL and ES solutions.
This degeneracy is similar to the “planet/binary-source” degeneracy, that is
often confronted in interpreting a short-term anomaly superposed on a single-
mass light curve, e.g., MOA-2012-BLG-486 (Hwang et al., 2013), in the sense
that explaining the observed anomaly requires one to include an additional
lens or source component. The third degeneracy type is caused by the ambiguity
between the cusp-approach solution and cusp-crossing solution, i.e., the
degeneracy between ESESin and ESESout solutions. This degeneracy arises
because the exact trajectory of the tertiary source ($S_{3}$) cannot be
specified because the region of the minor anomaly is not densely covered. A
similar degeneracy was reported in the interpretation of the short anomaly
appeared in the OGLE-2015-BLG-1459 light curve (Hwang et al., 2018). Finally,
the ambiguity of $S_{1}$ trajectory relative to the planetary caustic causes
an additional degeneracy in the interpretation of the minor anomaly. i.e,. the
degeneracy between the ESELin and ESELout solutions. This degeneracy was
mentioned by Gaudi & Gould (1997) for a single-host planetary event. The
inspection of the degeneracy types reveals that the origins of all the
degeneracy types arising in the interpretation of the event were already known
before. This implies that checking known types of degeneracies in analyzing
lensing light curves, especially with complex features, is important to
correctly characterize the lens system.

Figure 10.—  Middle panel: cumulative distributions of
$\Delta\chi^{2}=\chi^{2}_{\rm 2L1S}-\chi^{2}$ for the five degenerate models
from the base 2L1S model. The light curve in the top panel is inserted to show
where the fit improves. The two dotted lines indicate the times of the major
($t_{1}$) and minor ($t_{2}$) anomalies. Bottom panel: cumulative
distributions of the $\chi^{2}$ differences from the best-fit solution, i.e.,
$\Delta\chi^{2}=\chi^{2}-\chi^{2}_{\rm ESEL_{out}}$.

##  6\. Einstein radius

Measurement of $\theta_{\rm E}$ requires one to estimate the source radius,
i.e., $\theta_{\rm E}=\theta_{*}/\rho$. We derive $\theta_{*}$ from the source
color and brightness measured by regressing the observed $V$ and $I$
photometric data with the variation of the event magnification. We then
estimate de-reddened values, $(V-I,I)_{0}$, by calibrating the instrumental
values, $(V-I,I)$. This calibration is done utilizing the Yoo et al. (2004)
method, in which the reference position in the color-magnitude diagram (CMD)
for the calibration is the centroid of red giant clump (RGC) with
$(V-I,I)_{{\rm RGC},0}=(1.060,14.347)$ known by Bensby et al. (2013) and Nataf
et al. (2013).

In Figure 11, we mark the source (black dot) in the CMD. We also mark the RGC
centroid and blend. We will discuss in more detail about the nature of the
blend in Section 7. We note that the source position is determined using the
ELES solution for the reason to be discussed in Sect 7. The source positions
based on the other solutions are marked by grey-tone points. The measured
color and magnitude are $(V-I,I)=(1.703\pm 0.024,20.963\pm 0.003)$ for the
source and $(V-I,I)_{\rm RGC}=(1.978,15.688)$ for the RGC centroid. With the
measured offsets in the source color, $\Delta(V-I)$, and magnitude, $\Delta
I$, from the RGC centroid, the color and magnitude are calibrated as

$\displaystyle{(V-I,I)_{0}}$ | $\displaystyle{{}=(V-I,I)_{{\rm RGC},0}+\Delta(V-I,I)}$
---|---
| $\displaystyle{{}=(0.786\pm 0.024,19.622\pm 0.003).}$
(1)

The measured values points out that the source is a main sequence of a late G
type. We note that the source is well below the brightness limit of Gaia
observation.

In order to estimate the source radius, we convert $V-I$ into $V-K$, and
interpolate $\theta_{*}$ from the $(V-K)$–$\theta_{*}$ relation. Here we use
the Bessell & Brett (1988) relation for the color conversion, and the Kervella
et al. (2004) relation to derive $\theta_{*}$. The source is estimated to have
a radius of

$\theta_{*}=0.41\pm 0.03~{}\mu{\rm as}.$ (2)

For $\theta_{*}$ estimation, the source is assumed to lie at $D_{\rm S}=d_{\rm
GC}/[\cos l+\sin l(\cos\phi/\sin\phi)]\sim 7.85$ kpc, where we adopt a
Galactocentric distance of $d_{\rm GC}=8.16$ kpc and a bulge bar orientation
angle of $\phi=40^{\circ}$. The error of the $\theta_{*}$ measurement is
estimated based on the error of the measured source color and adding 7% error
in quadrature to account for combined uncertainty resulting from the RGC
centroiding and the color-$\theta_{*}$ conversion. We note that the source
companion contributes little flux, $q_{F}/(1+q_{F})\sim 4\%$, to the combined
source flux, and thus its effect on the estimated $\theta_{*}$ is minimal.
Together with $\rho$ and $t_{\rm E}$, the measured $\theta_{*}$ yields the
Einstein radius and the relative lens-source proper motion of

$\theta_{\rm E}={\theta*\over\rho}=0.91\pm 0.07~{}{\rm mas},$ (3)

and

$\mu={\theta_{\rm E}\over t_{\rm E}}=7.59\pm 0.59~{}{\rm mas}~{}{\rm
yr}^{-1},$ (4)

respectively. The estimated value of $\theta_{\rm E}$ is substantially bigger
than typical value of $\sim 0.5$ mas for events with low-mass stellar lenses
lying about halfway between Earth and the bulge. This suggests that the lens
lies at a close distance.

Figure 11.—  Source, RGC centroid, and blend positions in the CMD. The source
position is determined based on the ELES solution and those based on the other
solutions are marked by grey-tone points.

##  7\. Resolving Degeneracy

We check whether the degeneracies among the solutions can be resolved. For
this, we compare the flux ratio between $S_{1}$ and $S_{2}$ measured from
modeling, $q_{F,{\rm model}}=F_{S,2}/F_{S,1}$, with the flux ratio predicted
from the radius ratio $\rho_{2}/\rho_{1}$, $q_{F,{\rm pred}}$. The radius
ratio $\rho_{2}/\rho_{1}$ is measurable for the solutions in which the major
anomaly is explained by the caustic crossing of $S_{2}$, i.e., ESXX models. We
note that the major anomaly for the ELES model is explained by the caustic
crossing of $S_{1}$ instead of $S_{2}$, and thus this method cannot be
applied. In order to estimate $q_{F,{\rm pred}}$ from the $\rho_{2}/\rho_{1}$
ratio, we first estimate the physical radius of $S_{2}$ by

$R_{*,S_{2}}=\left({\rho_{2}\over\rho_{1}}\right)R_{*,S_{1}}.$ (5)

The ESXX solutions result in similar values of the color and brightness, as
shown in Figure 11, and thus we use a common physical source radius of
$R_{*,S_{1}}\sim 0.85~{}R_{\odot}$, which is deduced from $(V-I)_{0}$ and
$I_{0}$. We then estimate the absolute $I$-band magnitudes of the source
stars, $M_{I,1}$ and $M_{I,2}$, corresponding to $R_{*,S_{1}}$ and
$R_{*,S_{2}}$ from the tables of Pecaut et al. (2012) and Pecaut & Mamajek
(2013), and compute the flux ratio between $S_{1}$ and $S_{2}$ as

$q_{F,{\rm pred}}=10^{-0.4(M_{I,2}-M_{I,1})}.$ (6)

We note that a similar method could, in principle, be applied to the flux
ratio $q_{F,3}=F_{S,3}/F_{S,1}$ for the models with three source stars, i.e.,
ESES solutions. However, the normalized radius of the tertiary source,
$\rho_{3}$, is measured for none of the ESES model, and thus the method cannot
be implemented.

Figure 12.—  Plot of $q_{F,{\rm pred}}$ versus $q_{F,{\rm model}}$ for the
models explaining the major anomaly by introducing an extra source. Color
codings represent points with $\leq 1\sigma$ (red), $\leq 2\sigma$ (yellow),
$\leq 3\sigma$ (green), $\leq 4\sigma$ (cyan), and $\leq 5\sigma$ (blue). The
oblique solid line in each panel represents the positions at which $q_{F,{\rm
pred}}=q_{F,{\rm model}}$.

Figure 12 shows the scatter plot of points in the MCMC chains on the
$q_{F,{\rm model}}$–$q_{F,{\rm pred}}$ parameter plane for the four tested
ESXX models. In the plot, the color coding represents points with $\leq
1\sigma$ (red), $\leq 2\sigma$ (yellow), $\leq 3\sigma$ (green), $\leq
4\sigma$ (cyan), and $\leq 5\sigma$ (blue). The oblique solid line in each
panel represents the positions at which $q_{F,{\rm pred}}=q_{F,{\rm model}}$.
The plots show that the hypothesis of the second source’s caustic crossing for
the origin of the major anomaly is rejected at more than $3\sigma$ level,
suggesting that the ESXX models are unlikely to be correct interpretations of
the event. With these solutions rejected, the ELES solutions remain as the
only viable interpretation of the event.

##  8\. Physical lens parameters

Recognizing that only the ELES models provide plausible interpretations of the
observed lensing data, the lens parameters of the lens mass, $M$, and
distance, $D_{\rm L}$, are estimated based on the observables of the ELES
solutions. In order to unambiguously determine these parameters, one should
measure both $\theta_{\rm E}$ and $\pi_{\rm E}$, where $\pi_{\rm E}$
represents the microlens parallax (hereafter parallax), i.e.,

$M={\theta_{\rm E}\over\kappa\pi_{\rm E}};\qquad D_{\rm L}={{\rm
au}\over\pi_{\rm E}\theta_{\rm E}+\pi_{\rm S}}.$ (7)

Here $\kappa=4G/(c^{2}{\rm au})$, $D_{\rm S}$ is the source distance, and
$\pi_{\rm S}={\rm au}/D_{\rm S}$. The Einstein radius is firmly measured,
i.e., Equation (3). An additional modeling conducted considering the parallax
effect indicates that secure determinations of the parallax parameters are
difficult. Although $\pi_{\rm E}$ is not measured, the lensing observables of
$t_{\rm E}$ and $\theta_{\rm E}$ are related as

$t_{\rm E}={\theta_{\rm E}\over\mu};\qquad\theta_{\rm E}=(\kappa M\pi_{\rm
rel})^{1/2},$ (8)

where $\pi_{\rm rel}={\rm au}(D_{\rm L}^{-1}-D_{\rm S}^{-1})$. Then, the lens
parameters can be estimated by conducting a Bayesian analysis with the prior
Galactic model, defining the mass function, physical, and dynamical
distributions.

Figure 13.—  Bayesian posterior distributions of $M_{1}$ and $D_{\rm L}$. The
red and blue curves represent the distributions corresponding to the close and
wide models, respectively.

The Bayesian analysis is carried out in two steps. The first step is producing
events from a Monte Carlo simulation using a Galactic model. The models used
for the analysis include Han & Gould (2003), Han & Gould (1995), Zhang et al.
(2020) models for the physical, dynamical distributions, and mass function,
respectively. See section 5 of Han et al. (2020b) for more details about the
models. In the second step, the posterior distributions of $M$ and $D_{\rm L}$
are constructed for the events with observables, i.e., $t_{\rm E}$ and
$\theta_{\rm E}$, in the 1$\sigma$ ranges of the estimated observables among
the events produced from the simulation.

The posterior distributions $M_{1}$ and $D_{\rm L}$ derived from the ELES
solutions are shown in Figure 13, where the red curve is the distribution
resulting from the close solution, while the blue curve is the one from the
wide solution. In Table 7, we summarize the estimated parameters of $M_{1}$,
$M_{2}$, $M_{3}$, $D_{\rm L}$, $d_{\perp,2}$, and $d_{\perp,3}$, where the
last two quantities are the separations of $M_{2}$ and $M_{3}$ from $M_{1}$.
The listed values are the medians of the posterior distributions with errors
determined as 16% and 84% ranges. According to the estimated parameters from
the close model, the lens is a planetary system, in which a giant planet with
a mass

$M_{p}=2.56^{+1.13}_{-1.16}~{}M_{\rm J}$ (9)

belongs to a stellar binary composed of K and M dwarfs with masses

$M_{1}=0.61^{+0.27}_{-0.28}~{}M_{\odot},$ (10)

and

$M_{2}=0.15^{+0.07}_{-0.07}~{}M_{\odot},$ (11)

respectively. As predicted by the large $\theta_{\rm E}$, the lens is located
at a close distance of

$D_{\rm L}=3.86^{+1.00}_{-1.12}~{}{\rm kpc},$ (12)

suggesting that it is a disk object. The lens parameters derived from the wide
model are also listed in Table 7. According to the ELES models, the lens is
the seventh microlensing system with a planet in a stellar binary, followed by
OGLE-2006-BLG-284L (Bennett et al., 2020), OGLE-2007-BLG-349L (Bennett et al.,
2016), OGLE-2008-BLG-092L (Poleski et al., 2014), OGLE-2013-BLG-0341L (Gould
et al., 2014), OGLE-2016-BLG-0613L (Han et al., 2017), and OGLE-2018-BLG-1700L
(Han et al., 2020a).

Table 7Physical lens parameters Parameter | Close | Wide
---|---|---
$M_{1}$ ($M_{\odot}$) | $0.61^{+0.27}_{-0.28}$ | $0.59^{+0.26}_{-0.26}$
$M_{2}$ ($M_{\odot}$) | $0.15^{+0.07}_{-0.07}$ | $0.27^{+0.12}_{-0.12}$
$M_{3}$ ($M_{\rm J}$) | $2.56^{+1.13}_{-1.16}$ | $1.43^{+0.62}_{-0.62}$
$D_{\rm L}$ (kpc) | $3.86^{+1.00}_{-1.12}$ | $3.77^{+0.98}_{-1.13}$
$d_{\perp,2}$ (au) | $1.73^{+0.43}_{-0.48}$ | $7.82^{+2.04}_{-2.35}$
$d_{\perp,3}$ (au) | $3.32^{+0.82}_{-0.92}$ | $2.80^{+0.73}_{-0.84}$

The estimated physical lens parameters suggest the possibility that the
majority of the blended light comes from the lens. In Figure 11, We mark the
blend, with $(V-I,I)_{b}\sim(1.6,18.8)$, in the CMD. According to the
estimated $M$ and $D_{\rm L}$, the predicted brightness and color of the
primary lens, which explains most of the flux from the lens, are in the ranges
of

$18.3\lesssim I_{\rm L}\lesssim 22.3$ (13)

and

$1.3\lesssim(V-I)_{\rm L}\lesssim 2.9.$ (14)

Here we compute the lens brightness and color by $I_{\rm L}=M_{I}+5\log D_{\rm
L}-5+A_{I}$ and $(V-I)_{\rm L}=(V-I)_{0}+E(V-I)$, where $M_{I}$ and
$(V-I)_{0}$ are the absolute $I$-band magnitude and intrinsic color
corresponding to $M_{1}$, and we assume $A_{I}\sim 0.5$ and $E(V-I)\sim 0.4$
by adopting the half of the values to the source considering $D_{\rm L}$.
Then, the color and brightness of the blend are within the predicted ranges of
the lens.

We further check this hypothesis by measuring the offset between the centroid
of the apparent source image at the baseline and the source position in the
difference image during the lensing magnification. We measure

$(\Delta x,\Delta y)=(0.21\pm 0.10,0.19\pm 0.18)~{}{\rm arcsec}.$ (15)

Considering that the measured offset is within $\sim 2\sigma$ of the
measurement error, the hypothesis cannot be ruled out. Therefore, an important
portion of the blended light may come from the lens, but this can only be
established using high-resolution images that can be obtained from future
observations with the use of adaptive optics (AO) instrument mounted on
8m-class telescopes or space-based telescopes.

##  9\. Summary

We present the result from the investigation of KMT-2019-BLG-1715. The event
light curve displayed two anomalies from a typical caustic-crossing binary-
lensing light curve. We identified five pairs of solutions, in which the
anomalies were explained with the inclusion of an extra lens or source
component in addition to the base binary-lens model. We presented detailed
analysis for the individual solutions, and trace the origins of degeneracies.
To resolve the degeneracies, we compare the measured $S_{1}/S_{2}$ flux ratio
with the ratio deduced from the ratio of the source radii. Applying this
method left only a single pair of viable solutions, in which the anomaly with
a large deviation was produced by a third body of the lens, and the anomaly
with a small deviation was generated by a second source. A Bayesian analysis
indicated that the lens comprised a $\sim 2.6~{}M_{\rm J}$ planet and binary
stars with K and M spectral types lying in the disk. We pointed out that the
lens might be the blend, and this hypothesis could be confirmed from future
high-resolution followup observations.

Work by C.H. was supported by the grants of National Research Foundation of
Korea (2020R1A4A2002885 and 2019R1A2C2085965) Work by A.G. was supported by
JPL grant 1500811. The OGLE project has received funding from the National
Science Centre, Poland, grant MAESTRO 2014/14/A/ST9/00121 to A.U.. The MOA
project was supported by the grant 19KK0082 and 20H04754. Work by Y.H, was
supported by JSPS KAKENHI grant No. 17J02146. D.P.B., A.B., and C.R. were
supported by NASA through grant NASA-80NSSC18K0274. Work by N.K. was supported
by JSPS KAKENHI grant No. JP18J008. Work by C.R. was supported by an
appointment to the NASA Postdoctoral Program at the Goddard Space Flight
Center, administered by USRA through a contract with NASA. This research has
made use of the KMTNet system operated by the Korea Astronomy and Space
Science Institute (KASI) and the data were obtained at three host sites of
CTIO in Chile, SAAO in South Africa, and SSO in Australia.

## References

  * Alard & Lupton (1998) Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325
  * Albrow (2017) Albrow, M. 2017, MichaelDAlbrow/pyDIA: Initial Release on Github,Version v1.0.0, Zenodo, doi:10.5281/zenodo.268049
  * Albrow et al. (2001) Albrow, M. D., An, J., Beaulieu, J.-P., et al. 2001, ApJ, 549, 759
  * Alcock et al (1997) Alcock, C., Allsman, R. A., Alves, D., et al. 1997, ApJ, 479, 119
  * An & Han (2002) An, J. H., & Han, C. 2002, ApJ, 573, 351
  * Bennett (2010) Bennett, D. P. 2010, ApJ, 716, 1408
  * Bennett et al. (2016) Bennett, D. P., Rhie, S. H., Udalski, A., et al. 2016, AJ, 152, 125
  * Bennett et al. (2020) Bennett, D. P., Udalski, A., Bond, I. A., et al. 2020, AJ, 160, 72
  * Bensby et al. (2013) Bensby, T., Yee, J. C., Feltzing, S., et al. 2013, A&A, 549, 147
  * Bessell & Brett (1988) Bessell, M. S., & Brett, J. M. 1988, PASP, 100, 1134
  * Bond et al. (2001) Bond, I. A., Abe, F., Dodd, R. J., et al. 2001, MNRAS, 327, 868
  * Bond et al. (2002) Bond, I. A., Rattenbury, N. J., Skuljan, J., et al. 2002, MNRAS, 333, 71
  * Bozza (1999) Bozza, V. 1999, A&A, 348, 311
  * Bozza et al. (2018) Bozza, V., Bachelet, E., Bartolić, F., Heintz, T. M., Hoag, A. R., & Hundertmark, M. 2018, MNRAS, 479, 5157
  * Bozza et al. (2012) Bozza, V., Dominik, M., Rattenbury, N. J., et al. 2012, MNRAS, 424, 902
  * Daněk & Heyrovský (2015) Daněk, K., & Heyrovský, D. 2015, ApJ, 806, 99
  * Daněk & Heyrovský (2019) Daněk, K., & Heyrovský, D. 2019, ApJ, 880, 72
  * Di Stefano & Mao (1996) Di Stefano, R., & Mao, S. 1996, ApJ, 457, 93
  * Dominik (1999) Dominik, M. 1999, A&A, 349, 108
  * Dong et al. (2009) Dong, S., Bond, I. A., Gould, A., et al. 2009, ApJ, 698, 1826
  * Gaudi (1998) Gaudi, B. S. 1998, ApJ, 506, 533
  * Gaudi et al. (1998) Gaudi, B. S., Naber, R. M., & Sackett, P. D. 1998, ApJ, 502, L33
  * Gaudi & Gould (1997) Gaudi, B. S., & Gould A. 1997, ApJ, 486, 85
  * Gaudi & Han (2004) Gaudi, B. S., & Han, C. 2004, ApJ, 611, 528
  * Gould (1992) Gould, A. 1992, ApJ, 392, 442
  * Gould & Gaucherel (1997) Gould, A., & Gaucherel, C. 1997, ApJ, 477, 580
  * Gould et al. (2014) Gould, A., Udalski, A., Shin, I.-G., et al. 2014, Science, 345, 46
  * Griest & Safizadeh (1998) Griest, K., & Safizadeh, N. 1998, ApJ, 500, 37
  * Han et al. (2001) Han, C., Chang, H.-Y., An, J. H., & Chang, K. 2001, MNRAS, 328, 986
  * Han & Gould (1995) Han, C., & Gould, A. 1995, ApJ, 447, 53
  * Han & Gould (2003) Han, C., & Gould, A. 2003, ApJ, 592, 172
  * Han et al. (2021) Han, C., Lee, C.-U., Ryu, Y.-H., et al. 2021, A&A, in press [arXiv:2102.01806]
  * Han et al. (2020a) Han, C., Lee, C.-U., Udalski, A., et al. 2020a, AJ, 159, 48
  * Han et al. (2020b) Han, C., Shin, I.-G., Jung, Y. K., et al. 2020b, A&A, 641, A105
  * Han et al. (2017) Han, C., Udalski, A., Gould, A., et al. 2017, AJ, 154, 223
  * Hwang et al. (2013) Hwang, K. -H., Choi, J. -Y., Bond, I. A., et al. 2013 ApJ, 778, 55
  * Hwang et al. (2018) Hwang, K.-H., Udalski, A., Bond, I. A., et al. 2018, AJ, 155, 259
  * Kervella et al. (2004) Kervella, P., Thévenin, F., Di Folco, E., & Ségransan, D. 2004, A&A, 426, 29
  * Kim et al. (2016) Kim, S.-L., Lee, C.-U., Park, B.-G., et al. 2016, JKAS, 49, 37
  * Mao & Paczyński (1991) Mao, S., & Paczyński, B. 1991, ApJ, 374, 37
  * Nataf et al. (2013) Nataf, D. M., Gould, A., Fouqué, P., et al. 2013, ApJ, 769, 88
  * Pecaut & Mamajek (2013) Pecaut, M. J. & Mamajek, E. E. 2013, ApJS, 208, 9
  * Pecaut et al. (2012) Pecaut, M. J., Mamajek, E. E., & Bubar, E. J. 2012, ApJ, 746, 154
  * Poleski et al. (2014) Poleski, R., Skowron, J., Udalski, A., et al. 2014, ApJ, 795, 42
  * Ryu et al. (2010) Ryu, Y. -H., Han, C., Hwang, K.-H., et al. 2010, ApJ, 723, 81
  * Udalski (2003) Udalski, A. 2003, Acta Astron., 53, 291
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1
  * Udalski et al. (1994a) Udalski, A., Szymański, M., Kałużny, J., Kubiak, M., Mateo, M., & Krzemiński, W. 1994, Acta Astron. 44, 1
  * Udalski et al. (1994b) Udalski, A., Szymański, M., Mao, S., Di Stefano, R., Kałużny, J., Kubiak, M., Mateo, M., & Krzemiński, W. 1994, ApJ, 436, L103
  * Woźniak (2000) Woźniak, P. R. 2000, Acta Astron., 50, 42
  * Yee et al. (2012) Yee, J. C., Shvartzvald, Y., Gal-Yam, A., et al. 2012, ApJ, 755, 102
  * Yoo et al. (2004) Yoo, J., DePoy, D. L., Gal-Yam, A., et al. 2004, ApJ, 603, 139
  * Zhang et al. (2020) Zhang, X., Zang, W., Udalski, A., et al. 2020, AJ, 159, 116