# Chiral-Symmetric Higher-Order Topological Phases of Matter:
Supplemental Material

Wladimir A. Benalcazar<EMAIL_ADDRESS>Department of Physics, Princeton
University, Princeton, New Jersey 08542, USA  Department of Physics,
Pennsylvania State University, University Park, Pennsylvania 16802, USA
Alexander Cerjan  Sandia National Laboratories, Albuquerque, New Mexico 87185,
USA  Center for Integrated Nanotechnologies, Sandia National Laboratories,
Albuquerque, New Mexico 87185, USA

This Supplemental Material is divided into 4 sections. In Section I, we
describe the structure of chiral-symmetric Hamiltonians and how to calculate
the _winding number_ in discrete periodic extended systems (i.e., lattice
systems) in terms of a sublattice dipole moment operator, which will allow us
to derive the general class of multipole chiral numbers in any dimension. In
Section II, we enunciate the generalizations of these invariants for 2D and
3D, which label topological phases with higher-order topology. We also show
that these invariants are strictly quantized under chiral symmetry, and
provide expressions for them in the form of Bott index. Finally, we
differentiate these invariants from other first-order topological indices that
protect zero-dimensional states in 2D crystals. In Section III, we detail
properties of the model proposed in the Main Text, and show the correspondence
between the index $N_{xy}$ and its corner states. Finally, in Section IV, we
comment on the existence of bulk and boundary obstructions separating these
phases.

##  I Chiral-symmetric systems: general characteristics

Chiral-symmetric systems are described by Hamiltonians $\mathcal{H}$ that obey

$\displaystyle\Pi\mathcal{H}\Pi=-\mathcal{H},$ (S1)

where $\Pi$ is the _chiral operator_. The degrees of freedom in chiral-
symmetric systems can be divided into two sublattices, $A$ and $B$. The chiral
operator is equal to

$\displaystyle\Pi=\sum_{{\bf R},\alpha\in A}c^{\dagger}_{{\bf
R},\alpha}\left|0\right>\left<0\right|c_{{\bf R},\alpha}-\sum_{{\bf
R},\beta\in B}c^{\dagger}_{{\bf R},\beta}\left|0\right>\left<0\right|c_{{\bf
R},\beta}.$ (S2)

From now on, we will constrain our problem to the case in which the two
sublattices have equal number of degrees of freedom, $N_{A}=N_{B}$. For an
eigenstate of the Hamiltonian $\left|\psi_{n}\right>$ with energy
$\epsilon_{n}$, there is a partner eigenstate $\Pi\left|\psi_{n}\right>$ with
opposite energy,

$\displaystyle\mathcal{H}\Pi\left|\psi_{n}\right>=-\Pi\mathcal{H}\left|\psi_{n}\right>=-\epsilon_{n}\Pi\left|\psi_{n}\right>.$
(S3)

Thus, the energy spectrum of a chiral-symmetric system is symmetric with
respect to zero energy.

In the basis in which the degrees of freedom are ordered so that those of
sublattice $A$ come first and then those of sublattice $B$, the chiral
operator takes the form $\Pi=\sigma_{z}\otimes I_{N_{A}\times N_{A}}$, where
$\sigma_{z}$ the third Pauli matrix and $I_{n\times n}$ is the identity matrix
of size $n$. In that same basis, the Hamiltonian takes the form

$\displaystyle\mathcal{H}=\begin{pmatrix}0&h\\\ h^{\dagger}&0\end{pmatrix}.$
(S4)

The eigenstates of $\mathcal{H}$ can be written as
$\left|\psi\right>=\frac{1}{\sqrt{2}}(\psi^{A}_{n},\psi^{B}_{n})^{T}$, where
$\psi_{A}$ and $\psi_{B}$ are normalized vectors that exist only in the $A$,
$B$ subspaces, respectively. The chiral partner state with opposite energy is
$\Pi\left|\psi_{n}\right>=\frac{1}{\sqrt{2}}(\psi^{A}_{n},-\psi^{B}_{n})^{T}$.
Evaluating
$\mathcal{H}^{2}\left|\psi_{n}\right>=\epsilon_{n}^{2}\left|\psi_{n}\right>$
leads to the eigenvalue problems

$\displaystyle(hh^{\dagger})\psi^{A}_{n}$
$\displaystyle=\epsilon_{n}^{2}\psi^{A}_{n}$
$\displaystyle(h^{\dagger}h)\psi^{B}_{n}$
$\displaystyle=\epsilon_{n}^{2}\psi^{B}_{n},$ (S5)

so that $\psi^{A}_{n}$ and $\psi^{B}_{n}$ can be easily obtained by
diagonalizing $hh^{\dagger}$ or $h^{\dagger}h$, respectively. Equation (S5)
has a structure compatible with a singular value decomposition. Indeed, it is
possible to write $h$ as

$\displaystyle h=U_{A}\Sigma U_{B}^{\dagger}$ (S6)

where $U_{A}$ is a unitary matrix representing the space spanned by
$\\{\psi^{A}_{n}\\}$, i.e.,
$U_{A}=(\psi^{A}_{1},\psi^{A}_{2}\ldots,\psi^{A}_{N_{A}})$ and similarly for
$U_{B}$, and $\Sigma$ is a diagonal matrix containing the singular values.
Notice, for example, that

$\displaystyle hh^{\dagger}=U_{A}\Sigma^{2}U_{A}^{\dagger}$ $\displaystyle
h^{\dagger}h=U_{B}\Sigma^{2}U_{B}^{\dagger}$ (S7)

which are compatible with Eq. (S5) if we identify the squared energies
$\\{\epsilon^{2}_{n}\\}$ with the square of the singular values in $\Sigma$.

In what follows, it will be convenient to have a ‘flattened’ version of
Hamiltonian (S4), given by

$\displaystyle\mathcal{H}=\begin{pmatrix}0&q\\\
q^{\dagger}&0\end{pmatrix},\quad q=U_{A}U^{\dagger}_{B},$ (S8)

whose energies are $\epsilon_{n}=\pm 1$. The flattening was achieved by
dropping the singular value matrix $\Sigma$ in the definition of $q$ in Eq.
(S8). Contrast this with the equivalent expression in Eq. (S6) for Hamiltonian
(S4). Although the Hamiltonians (S4) and (S8) represent different systems,
those two systems will belong to the same topological phase. We will see,
however, that Hamiltonian (S8) facilitates the calculation of the topological
invariants associated with its topological phase.

###  I.1 The winding number as a sublattice polarization

We now turn our attention to chiral-symmetric Hamiltonians that are extended
and periodic (i.e., lattice systems). These systems have topological phases
labelled by an integer invariant, known as the _winding number_. Since
topological properties of a topological phase are preserved under adiabatic
deformations, i.e., those deformations that preserve the energy gap and
symmetry, let us consider, without loss of generality, a flattened version of
Hamiltonian (S4) at each crystal momentum $k$,

$\displaystyle\mathcal{H}(k)=\begin{pmatrix}0&q(k)\\\
q(k)^{\dagger}&0\end{pmatrix},\quad q(k)=U_{A}(k)U^{\dagger}_{B}(k),$ (S9)

where $\mathcal{H}(k)$ is a Bloch Hamiltonian with a _flat_ energy spectrum
$\epsilon=\pm 1$ at all $k$. The winding number, which classifies different
homotopy classes $\pi_{1}[U(n)]\cong\mathbb{Z}$ in 1D, is given by

$\displaystyle
N_{x}=\frac{\mathrm{i}}{2\pi}\int_{-\pi}^{\pi}\text{Tr}\left[q(k)^{\dagger}\partial_{k}q(k)\right]dk\in\mathbb{Z}.$
(S10)

By replacing $q(k)=U_{A}(k)U^{\dagger}_{B}(k)$ in Eq. (S10), we have

$\displaystyle N_{x}$
$\displaystyle=\frac{\mathrm{i}}{2\pi}\int_{-\pi}^{\pi}\text{Tr}\left[U_{A}^{\dagger}(k)\partial_{k}U_{A}(k)\right]dk+\frac{\mathrm{i}}{2\pi}\int_{-\pi}^{\pi}\text{Tr}\left[U_{B}(k)\partial_{k}U_{B}^{\dagger}(k)\right]dk.$
$\displaystyle=-\frac{1}{2\pi}\int_{-\pi}^{\pi}\text{Tr}\left[\mathcal{A}_{A}(k)\right]dk+\frac{1}{2\pi}\int_{-\pi}^{\pi}\text{Tr}\left[\mathcal{A}_{B}(k)\right]dk.$
(S11)

Here $\mathcal{A}_{A}(k)=-\mathrm{i}U_{A}^{\dagger}(k)\partial_{k}U_{A}(k)$ is
the Berry connection in the $A$ sublattice, as is $\mathcal{A}_{B}$ for the
$B$ sublattice. In the second step we have used
$\text{Tr}[U^{\dagger}(k)\partial_{k}U(k)]=-\text{Tr}[U(k)\partial_{k}U^{\dagger}(k)]$.
Notice that Eq. (S11) is _not_ defined mod 1.

###  I.2 Winding number in discrete space

We now seek to derive an expression equivalent to the winding number but in
real space instead of momentum space. Such an expression would enable the
characterization of disordered lattices. A first step is to discretize Eq.
(S11). For that purpose, consider the following expression,

$\displaystyle U^{\dagger}(k)\partial_{k}U(k)\approx U^{\dagger}(k)(U(k+\delta
k)-U(k))/\delta k,$

valid as $\delta k\rightarrow 0$. Rearranging this expression,

$\displaystyle I+\delta kU^{\dagger}(k)\partial_{k}U(k)\approx
U^{\dagger}(k)U(k+\delta k).$

Taking the logarithm on both sides, and using the approximation
$\text{Log}(1+x)\approx x$ on the LHS, we have

$\displaystyle\delta
kU^{\dagger}(k)\partial_{k}U(k)\approx\text{Log}(U^{\dagger}(k)U(k+\delta
k)).$ (S12)

Let us now use this expression on Eq. (S11) to have a discretized version of
it,

$\displaystyle N_{x}$
$\displaystyle=\frac{\mathrm{i}}{2\pi}\sum_{k}\text{Tr}\left[\text{Log}(U_{A,k}^{\dagger}U_{A,k+\delta
k})-\text{Log}(U_{B,k}^{\dagger}U_{B,k+\delta k})\right]$
$\displaystyle=\frac{1}{2\pi\mathrm{i}}\sum_{k}\text{Tr}\text{Log}\left[F_{A,k}F_{B,k}^{\dagger}\right],$
(S13)

where in the last step we used the fact that, for a unitary matrix $M$,
$\text{Log}(M)=-\text{Log}(M^{\dagger})$. The unitary matrices
$F_{\mathcal{S},k}$ are determined by first considering the matrices

$\displaystyle G_{\mathcal{S},k}=U_{\mathcal{S},k+\delta
k}^{\dagger}U_{\mathcal{S},k},$ (S14)

which are _not_ unitary due to the discretization of $k$ 111$G_{k}$ becomes
unitary only in the thermodynamic limit, where the spectrum $k$ is
continuous.. To restore unitatiry, consider the singular value decomposition

$\displaystyle
G_{\mathcal{S},k}=V_{\mathcal{S},k}D_{\mathcal{S},k}W_{\mathcal{S},k}^{\dagger},$
(S15)

where $D_{\mathcal{S},k}$ is the diagonal matrix containing the singular
values. In the thermodynamic limit, $D_{\mathcal{S},k}$ is the identity
matrix. For finite $L$, on the other hand, the magnitud of the diagonal
elements of $D_{\mathcal{S},k}$ are less than 1. We define

$\displaystyle F_{\mathcal{S},k}=V_{\mathcal{S},k}W_{\mathcal{S},k}^{\dagger}$
(S16)

which is unitary.

###  I.3 Winding number in position space

Having discretized the expression for the winding number for the case in which
the momenta is discrete, we now derive the position space version of the
winding number [2]. For that purpose, consider the second-quantization
operators for the anihilation of electrons in real space and crystal momentum
space, which are related by the Fourier transform

$\displaystyle c_{R,\alpha}$
$\displaystyle=\frac{1}{\sqrt{L}}\sum_{k}e^{-\mathrm{i}kR}c_{k,\alpha},$
$\displaystyle c_{k,\alpha}$
$\displaystyle=\frac{1}{\sqrt{L}}\sum_{R}e^{\mathrm{i}kR}c_{R,\alpha}.$ (S17)

To enforce periodic boundary conditions, we impose the constraint

$\displaystyle c_{R+L,\alpha}=c_{R,\alpha}\rightarrow
k=\frac{2\pi}{L}m,\;\;m\in\mathbb{Z},$ (S18)

which means that in Eq. (S17), $k\in\delta_{k}(0,1,\ldots L-1)$, where
$\delta_{k}=2\pi/L$.

To determine the winding number, we will need to make use of the operator
associated with polarization [3] defined over each sublattice
$\mathcal{S}=A,B$ [2],

$\displaystyle\hat{P}_{x}^{\mathcal{S}}=\sum_{R,\alpha\in\mathcal{S}}c^{\dagger}_{R,\alpha}\left|0\right>\exp\left[-\mathrm{i}\frac{2\pi}{L}R\right]\left<0\right|c_{R,\alpha}.$
(S19)

This operator can be written in momentum space as

$\displaystyle\hat{P}_{x}^{\mathcal{S}}=\sum_{k,\alpha\in\mathcal{S}}c^{\dagger}_{k+\delta_{k},\alpha}\left|0\right>\left<0\right|c_{k,\alpha}.$
(S20)

The action of $\hat{P}_{x}^{\mathcal{S}}$ on a state
$\left|k,\alpha\right>=c^{\dagger}_{k,\alpha}\left|0\right>$ is to shift the
momentum of the state on sublattice $\mathcal{S}$, i.e.,
$\hat{P}_{x}^{\mathcal{S}}\left|k,\alpha\right>=\left|k+\delta_{k},\alpha\right>$.

In what follows, we show that the most general expression for the winding
number is

$\displaystyle\boxed{N_{x}=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(\bar{P}_{x}^{A}\bar{P}_{x}^{B\dagger}\right)}$
(S21)

where
$\bar{P}_{x}^{\mathcal{S}}=U_{\mathcal{S}}^{\dagger}P_{x}^{\mathcal{S}}U_{\mathcal{S}}$.

For that purpose, we will transform Eq. (S21) into crystal momentum space and
show that it reduces to Eq. (S13). For that purpose, let us consider the
projectors into sublattice $\mathcal{S}$

$\displaystyle\mathcal{P}_{\mathcal{S}}=\sum_{n,k,\alpha,\beta}[u^{\mathcal{S}}_{n,k}]_{\alpha}c^{\dagger}_{k,\alpha}\left|0\right>\left<0\right|c_{k,\beta}[u^{\mathcal{S}}_{n,k}]^{*}_{\beta},$
(S22)

where $[u^{\mathcal{S}}_{n,k}]_{\alpha}$ denotes the $\alpha$th component of
the $n$th singular state in the SVD decomposition of Eq. (S6), corresponding
to the $n$th column of $U_{\mathcal{S}}$, for $\mathcal{S}=A$ or $B$, and ∗
denotes complex conjugation. Written compactly,
$\mathcal{P}_{\mathcal{S}}=U_{\mathcal{S}}U_{\mathcal{S}}^{\dagger}$.
Similarly, we can spell out $q=U_{A}U_{B}^{\dagger}$ [Eq. (S8)] as

$\displaystyle
q=\sum_{n,k,\alpha,\beta}[u^{A}_{n,k}]_{\alpha}c^{\dagger}_{k,\alpha}\left|0\right>\left<0\right|c_{k,\beta}.[u^{B}_{n,k}]^{*}_{\beta}$
(S23)

We need to calculate
$\mathcal{P}_{A}P^{A}_{x}qP^{B\dagger}_{x}\mathcal{P}_{B}=U_{A}U_{A}^{\dagger}P^{A}_{x}U_{A}U_{B}^{\dagger}P^{B\dagger}_{x}U_{B}U_{B}^{\dagger}$.
First, let us calculate

$\displaystyle\mathcal{P}_{A}P_{x}^{A}=\sum_{n,k,\alpha}\left|u^{A}_{n,k+\delta
k}\right>[u^{A}_{n,k+\delta k}]^{*}_{\alpha}\left<0\right|c_{k,\alpha}$
$\displaystyle
P_{x}^{B\dagger}\mathcal{P}_{B}=\sum_{n,k,\alpha}c^{\dagger}_{k,\alpha}\left|0\right>[u^{B}_{n,k+\delta
k}]^{*}_{\alpha}\left<u^{B}_{n,k+\delta k}\right|,$ (S24)

where
$\left|u^{\mathcal{S}}_{n,k}\right>=\sum_{\alpha}[u^{\mathcal{S}}_{n,k}]_{\alpha}c^{\dagger}_{k,\alpha}\left|0\right>$.
Using Eqs. S23 and S24, we have

$\displaystyle\mathcal{P}_{A}P^{A}_{x}qP_{x}^{B\dagger}\mathcal{P}_{B}$
$\displaystyle=\sum_{k,n,m,l}\sum_{\alpha,\beta}\left|u^{A}_{n,k+\delta
k}\right>[u^{A}_{n,k+\delta
k}]^{*}_{\alpha}[u^{A}_{m,k}]_{\alpha}[u^{B}_{m,k}]^{*}_{\beta}[u^{B}_{l,k+\delta
k}]_{\beta}\left<u^{B}_{l,k+\delta k}\right|$
$\displaystyle=\sum_{k,n,m,l}\left|u^{A}_{n,k+\delta
k}\right>[F_{A,k}]_{n,m}[F^{\dagger}_{B,k}]_{m,l}\left<u^{B}_{l,k+\delta
k}\right|$ $\displaystyle=\sum_{k,n,l}\left|u^{A}_{n,k+\delta
k}\right>[F_{A,k}F^{\dagger}_{B,k}]_{n,l}\left<u^{B}_{l,k+\delta k}\right|$
$\displaystyle=\sum_{k,n,l}\left|u^{A}_{n,k}\right>[F_{A,k-\delta
k}F^{\dagger}_{B,k-\delta k}]_{n,l}\left<u^{B}_{l,k}\right|.$ (S25)

The expression (S25) implies that
$\bar{P}_{x}^{A}\bar{P}_{x}^{B\dagger}=U_{A}^{\dagger}P^{A}_{x}U_{A}U_{B}^{\dagger}P_{x}^{B\dagger}U_{B}$
takes the form

$\displaystyle\bar{P}_{x}^{A}\bar{P}_{x}^{B\dagger}=\left(\begin{array}[]{ccccc}F_{A,k_{0}}F^{\dagger}_{B,k_{0}}&0&\ldots&0\\\
0&F_{A,k_{1}}F^{\dagger}_{B,k_{1}}&\ldots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\ldots&F_{A,k_{L-1}}F^{\dagger}_{B,k_{L-1}}\\\ \end{array}\right).$ (S30)

Notice that the matrix is diagonal in $k$. Therefore, we have

$\displaystyle N_{x}$
$\displaystyle=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(\bar{P}_{x}^{A}\bar{P}_{x}^{B\dagger}\right)$
$\displaystyle=\frac{1}{2\pi\mathrm{i}}\sum_{k}\text{Tr}\text{Log}\left[F_{A,k}F_{B,k}^{\dagger}\right]$
(S31)

which is identical to the expression for the winding number (S13).

In contrast, the invariant that protects the quantization of the dipole moment
in chiral-symmetric insulators at half-filling is

$\displaystyle
p_{x}=-\frac{1}{2\pi}\text{Im}\text{Log}[\text{Det}(U_{\textrm{occ}}^{\dagger}P_{x}U_{\textrm{occ}})\text{Det}(P_{x}^{\dagger
1/2})],$ (S32)

where $U_{\textrm{occ}}$ is the subspace of occupied bands at half filling and
$P_{x}=\sum_{R,\alpha}\left|R,\alpha\right>\mathcal{U}_{D}(R)\left<R,\alpha\right|$
is the position operator defined over the entire lattice [3]. The dipole
moment invariant is quantized to $p_{x}=0$ or $1/2$ under chiral symmetry, and
is related to the winding number by the expression $p_{x}=N_{x}/2$ mod 1.

##  II Generalizations of the winding number to higher-order topological
systems with chiral symmetry

As stated in the Main Text, we generalize the invariant (S21) to 2D and 3D. In
2D, the invariant is

$\displaystyle\boxed{N_{xy}=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(\bar{Q}_{xy}^{A}\bar{Q}_{xy}^{B\dagger}\right)\in\mathbb{Z}}$
(S33)

where
$\bar{Q}_{xy}^{\mathcal{S}}=U_{\mathcal{S}}^{\dagger}Q^{\mathcal{S}}_{xy}U_{\mathcal{S}}$
for $\mathcal{S}=A,B$, and $Q_{xy}^{\mathcal{S}}$ is the sublattice quadrupole
moment operator

$\displaystyle Q^{\mathcal{S}}_{xy}$ $\displaystyle=\sum_{{\bf
R},\alpha\in\mathcal{S}}\left|{\bf
R},\alpha\right>\text{Exp}\left(-\mathrm{i}\frac{2\pi
xy}{L_{x}L_{y}}\right)\left<{\bf R},\alpha\right|.$ (S34)

Similarly, in 3D the invariant is

$\displaystyle\boxed{N_{xyz}=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(\bar{O}_{xyz}^{A}\bar{O}_{xyz}^{B\dagger}\right)\in\mathbb{Z}}$
(S35)

where
$\bar{O}^{\mathcal{S}}_{xyz}=U_{\mathcal{S}}^{\dagger}O^{\mathcal{S}}_{xyz}U_{\mathcal{S}}$,
for $\mathcal{S}=A,B$, and $O_{xyz}^{\mathcal{S}}$ is the sublattice octupole
moment operator

$\displaystyle O^{\mathcal{S}}_{xyz}$ $\displaystyle=\sum_{{\bf
R},\alpha\in\mathcal{S}}\left|{\bf
R},\alpha\right>\text{Exp}\left(-\mathrm{i}\frac{2\pi
xyz}{L_{x}L_{y}L_{z}}\right)\left<{\bf R},\alpha\right|,$ (S36)

###  II.1 Quantization of the real space invariants

The invariants (S21), (S33), and (S35) take the form

$\displaystyle
N=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(U_{A}^{\dagger}M_{A}U_{A}U_{B}^{\dagger}M^{\dagger}_{B}U_{B}\right)$
(S37)

where $M_{\mathcal{S}}$ (for $\mathcal{S}=A,B$) is $P_{x}^{\mathcal{S}}$,
$Q_{xy}^{\mathcal{S}}$, or $O_{xyz}^{\mathcal{S}}$ in 1D, 2D, or 3D,
respectively. Notice that since the matrices $M_{\mathcal{S}}$ and
$U_{\mathcal{S}}$ are unitary, we have

$\displaystyle\text{Det}(U_{A}^{\dagger}M_{A}U_{A}U_{B}^{\dagger}M_{B}^{\dagger}U_{B})=\text{Det}(M_{A}M^{\dagger}_{B})=1,$
(S38)

where in the last step we use the fact that $M_{x}^{A}=M_{x}^{B}$, which is
true because the two sublattices have (i) equal number of degrees of freedom
in each unit cell and (ii) the same number of unit cells, from which it
follows that $M_{A}$ and $M_{B}$ are identical. From (S38) it follows that
tracing the logarithm of
$U_{A}^{\dagger}M_{A}U_{A}U_{B}^{\dagger}M^{\dagger}_{B}U_{B}$ in (S37) will
necessarily have to give a phase that is a multiple of $2\pi i$, i.e., it will
be of the form $2\pi\mathrm{i}N$, with $N\in\mathbb{Z}$. This integer $N$ is
indeed the topological invariant.

###  II.2 Real space invariants written as a Bott index

Starting again with Eq. S37, and performing a unitary transformation using
$U_{B}$ of the matrix product on the interior of the matrix logarithm, we can
rewrite (S37) as

$N=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(U_{B}U_{A}^{\dagger}M_{A}U_{A}U_{B}^{\dagger}M^{\dagger}_{B}U_{B}U_{B}^{\dagger}\right).$
(S39)

Noting the definition of $q=U_{A}U_{B}^{\dagger}$ and again that
$M_{A}=M_{B}=M$ are unitary (see discussion in previous section),

$\displaystyle N$
$\displaystyle=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(q^{-1}M_{A}qM_{B}^{\dagger}\right)$
$\displaystyle=\frac{1}{2\pi\mathrm{i}}\text{Tr}\text{Log}\left(q^{-1}MqM^{-1}\right)$
$\displaystyle=\textrm{Bott}(q^{-1},M).$ (S40)

###  II.3 Topological protection of zero-dimensional states of first-order
topolgy

Chiral-symmetric systems can also protect zero-dimensional states at
topological defects. For a point defect in $d$ dimensions, $q({\bf k},{\bf
r})$ in Eq. (S4) is parametrized by $d$ momentum variables ${\bf k}$ and $d-1$
position variables ${\bf r}$, and the number of protected states is given by
the homotopy class $\pi_{2d-1}[U(N)]=\mathbb{Z}$ into which $q({\bf k},{\bf
r})$ falls [4], such as in the 2D Jackiw and Rossi lattice model for the
protection of 0D topological states bound to vortices. Nontrivial homotopies
exist only in odd dimensional manifolds, and thus rule out, for example, the
protection of 0D states at the corners of 2D crystals. In contrast, our
invariants, Eqs. (5) and (6) of the Main Text (Eq. (S33) and (S35) here),
allow for such protection, and thus are of different nature than first-order
topological indices based on homotopy classes.

##  III The extended QTI model

In this section we describe in detail certain characteristics of the QTI model
with long-range hopping terms (Eq. (7), (8), and (9) in the Main Text). We
show (A) how these phases are in general outside the framework of _topological
quantum chemistry_ [5, 6], (B) how the phase diagram is modified when the
horizontal hopping terms are different from the vertical ones, (C) the support
of the corner states, (D) the correspondence between the topological index
$N_{xy}$ [Eq. (5) in the Main Text, or Eq. (S33) here] and the number and
chiral charge of the corner states, (E) some subtleties regarding how the
phase diagrams in the Main Text were determined, and (F) the convergence of
our results with disorder with system size.

###  III.1 Symmetry representations of the topological phases

The extended QTI model has the Hamiltonian of Eq. (S4) with off-diagonal term

$\displaystyle h({\bf k})=h_{QTI}({\bf k})+C_{s}h_{\textrm{SLR}}({\bf
k})+C_{d}h_{\textrm{DLR}}({\bf k})$ (S41)

where $h_{\textrm{QTI}}({\bf k})$, $h_{\textrm{SLR}}({\bf k})$, and
$h_{\textrm{DLR}}({\bf k})$ are defined in Eqs. (7), (8), and (9) of the Main
Text.

Consider the case $C_{s}=1$, $C_{d}=0$. In that case, the Hamiltonian is
$C_{4v}$ symmetric and supports bulk-obstructed topological phases. At high-
symmetry points (HSPs) and lines of the BZ, the representations of the
elements of $C_{4v}$ are given in Table S1. To capture more directly the
topological protection due to crystalline symmetries, Table S1 also shows the
symmetry indicator invariants, defined as

$\displaystyle[M_{j}]=\\#M_{j}-\\#\Gamma_{j}$ (S42)

where $\\#M_{j}$ is the number of bands below the gap with $C_{4}$ symmetry
representations $M_{j}=e^{\mathrm{i}\pi(2j-1)/4}$, for $j=1,2,3,4$, at the
${\bf M}$ point of the BZ, and similarly for the ${\bf\Gamma}$ point.

phase | irreps at ${\bf\Gamma}$ | irreps at ${\bf M}$ | ($[M_{1}]$, $[M_{2}]$, $[M_{3}]$, $[M_{4}]$)
---|---|---|---
$N_{xy}=0$ | $\\{e^{\mathrm{i}3\pi/4},e^{-\mathrm{i}3\pi/4}\\}$ | $\\{e^{\mathrm{i}3\pi/4},e^{-\mathrm{i}3\pi/4}\\}$ | (0,0,0,0)
$N_{xy}=1$ | $\\{e^{\mathrm{i}3\pi/4},e^{-\mathrm{i}3\pi/4}\\}$ | $\\{e^{\mathrm{i}\pi/4},e^{-\mathrm{i}\pi/4}\\}$ | (1,-1,-1,1)
$N_{xy}=4$ | $\\{e^{\mathrm{i}3\pi/4},e^{-\mathrm{i}3\pi/4}\\}$ | $\\{e^{\mathrm{i}3\pi/4},e^{-\mathrm{i}3\pi/4}\\}$ | (0,0,0,0)
Table S1:  $C_{4}$ symmetry representations and symmetry indicator invariants
(S42) for the lowest two bands of Hamiltonian (S41) at the $C_{4}$ invariant
points of the BZ ${\bf\Gamma}$ and ${\bf M}$ for all topological phases. At
all $C_{2}$-invariant points, the representations are $\pm\mathrm{i}$.
Similarly, at all $M_{x}$ and $M_{y}$ invariant lines, the representations are
$\pm 1$. In both cases, those representations lead to trivial symmetry
indicator invariants for $C_{2}$ and reflection symmetries.

A result worth noting from Table S1 is that, while the $N_{xy}=1$ phase (i.e.,
the QTI nontrivial phase with $q_{xy}=e/2$) has nontrivial symmetry indicator
invariants, the $N_{xy}=4$ does not. Consequently, the topology of this last
phase cannot be captured by the theory of induction of representations and
symmetry indicator invariants [5, 6]. Nevertheless, its nontrivial nature is
evident in both the existence of a nonvanishing bulk invariant $N_{xy}=4$ and
the existence of an equal number of topological zero-energy states at each
corner [Fig. 1(d) of the Main Text and Fig. S3 here].

The representations due to $C_{2}$ symmetry and reflection symmetries are
equal at all HSPs of the BZ in all phases, and thus have trivial symmetry
indicator invariants.

###  III.2 Phase diagram for the ordered $C_{2v}$ symmetric system

In the interest of completeness, in Fig. S1 we show the phase diagram for a
system which is $C_{2v}$ symmetric, with $w_{m,x}=3w_{m,y}$, $w_{D}=0$, and
$w_{m\geq 0}=0$.

Figure S1: Phase diagram of the chiral-symmetric second-order topological
phase, $N_{xy}$ (left), and the quadrupole phase $q_{xy}$ (right), for a
$C_{2v}$ symmetric system with open boundaries as a function of the hopping
ratios, $w_{1,y}/v$ and $w_{2,y}/v$, with $w_{m,x}=3w_{m,y}$ and $w_{D}=0$.
Phase transitions where the boundary band gap closes along the $x$ edge are
shown in green, while those which close along the $y$ edge are shown in cyan.

###  III.3 Corner states and their support

In this section, we show the form and localization of two sets of corner
states for phases indicated in the Main Text. First, in Fig. S3, we show one
possible choice of the four corner-localized states of the lattice considered
in Fig. 1 in the Main Text with $N_{xy}=4$. As can be seen, the summation of
the probability density functions of the four orthogonal states equals that of
the total set of four states in a single corner, as shown in the Main Text.
Second, in Fig. S3 we show the probability density function of the corner
states of a $N_{xy}=-1$ phase from the system considered in Fig. 2 in the Main
Text. Notice that the corner states have support on the opposite sublattice as
those in the $N_{xy}=4$ phase (Fig. S3) due to the reversal in the sign of the
bulk index $N_{xy}$.

Figure S2: Probability density functions of the four corner states at one of
the corners of a $N_{xy}=4$ phase (four left plots) and cumulative PDF for all
the four states (rightmost plot). Support at each sublattice is indicated by
blue and red colors, respectively. All corner states in the same corner have
the same chiral charge and thus have support on only one sublattice. In this
simulation, $v=1$, $w_{1}=2$, $w_{2}=2$, $w_{D}=0$.

Figure S3: Probability density functions of all the corner states of a
$N_{xy}=-1$ phase (4 states in total, one per corner). Support at each
sublattice is indicated by blue and red colors, respectively. In this
simulation, $v=1$, $w_{1}=1$, $w_{2}=0.8$, $w_{D}=0.5$.

###  III.4 Correspondence between $N_{xy}$, band gap closings, and corner
states

To claim that the multipole chiral numbers introduced here are the topological
invariant connected with the appearance of corner-localized modes across a
band gap closing, one must show that all three phenomena are causally
connected. In other words, changes in either $N_{xy}$ or the number (or type)
of corner-localized modes necessarily imply a change in other quantity, and
moreover that these changes can only occur at a band gap closing.

Figure S4: Correspondence between the appearance of zero energy states, change
in the topological invariant $N_{xy}$, and the closing of the bulk and edge
band gaps for two of the systems considered in the Main Text in Figs. 1(b) and
2(a). (a) Copy of Fig. 1(b) from the Main Text, red line indicates the choice
of $w_{1}/v$ and $w_{2}/v$ in (b-d). (b) Number of states relative to their
energy and $w_{2}/v$ calculated with open boundary conditions. (c) $N_{xy}$ as
a function of $w_{2}/v$ calculated using open boundary conditions. (d) Bulk
band gap, $\Delta\epsilon$ as a function of $w_{2}/v$ calculated using
periodic boundary conditions. (e-h) Similar to (a-d), except for the system in
Fig. 2(a) of the Main Text. (i) Edge band gap calculated using a semi-open
boundary, i.e. periodic in $x$, but open in $y$. In all simulations of both
systems, the system size was $80\times 80$ unit cells.

To support this claim, Fig. S4(a-d) presents numerical simulations of all
three of these quantities, $N_{xy}$, the distribution of states in energy, and
the bulk band gap, for the long-range quadrupole topological insulator model
considered in Fig. 1(b) of the Main Text with fixed $w_{1}/v=0.4$. As can be
seen, the closing of the bulk band gap at $w_{2}/v=1$ coincides with 8 states
departing from each of the upper and lower bulk bands and becoming pinned at
$\epsilon=0$. These 16 states in total are those corner-localized states shown
in Fig. 1(d) of the Main Text, which are predicted to exist as $N_{xy}=4$ for
$w_{2}/v>1$. The slight discrepancy between the change in $N_{xy}$ and the
closing of the bulk band gap is due to finite system effects, where the
topological invariant does not change until the bulk band gap is larger that
the approximate coupling strength between the corner states in adjacent
corners.

In Fig. S4(e-i), we show this same correspondence for the quadrupole
topological insulator with both long-range straight and long-range diagonal
hoppings considered in Fig. 2(b) of the Main Text with fixed $w_{1}/v=0.8$. As
$w_{2}/v$ is increased from $0$ to $2$, this system undergoes three separate
phase transitions. The first transition, $N_{xy}=1\rightarrow-1$ near
$w_{2}/v=0.4$, occurs at a closing of the edge band gap, not the bulk band
gap. This particular phase transition also does not change the number of
corner-localized states, as there is only one state localized to each corner
in both phases. Rather, across the phase transition, the corner states swap
the sublattices on which they are supported. The other two phase transitions,
$N_{xy}=-1\rightarrow 0$ and $N_{xy}=0\rightarrow 4$, both occur when the bulk
band gap closes, and result in changes in the number of corner-localized
states. Again, discrepancies between the change in $N_{xy}$ from the band gap
closings are due to finite system size effects.

###  III.5 Determination of the phase diagrams of the extended QTI model

There are some subtleties in how the phase diagrams shown in Figs. 1(b) and
2(a) of the Main Text are constructed due to finite system size effects. As
was shown in Fig. S4, computing the topological invariants $N_{xy}$ in finite
systems causes the topological invariant not to change at exactly the same
location where the band gap closes. Moreover, it is not numerically feasible
to calculate an entire 2D phase diagram at a reasonable sampling density even
at the system sizes shown in Fig. S4 ($80\times 80$ unit cells), due to the
memory requirements necessary to calculate the full singular value
decomposition of $H$. In the definition of $N_{xy}$, every single singular
value is necessary to achieve accurate results.

As such, the phase diagrams shown in Figs. 1b, 2a, and 2b are constructed by
first finding closings in the bulk and edge band structures, assuming the
system is infinite in any non-open direction. Numerically, artifacts can
appear in the calculation of the band gap closings due to finite system sizes,
as shown in Fig. S5, but simulations strongly suggest that the phase diagrams
converge to those shown in the Main Text in the thermodynamic limit. Then, we
calculate $N_{xy}$ within each region at an ensemble of points using open
boundaries in all directions and using large system sizes. We also verify
these calculations through direct calculation of $N_{xy}$ across the entire 2D
parameter space but using smaller system sizes.

Figure S5: Bulk band gaps of the $C_{4v}$ symmetric system considered in Fig.
1(b) of the Main Text calculated using periodic boundary conditions with the
spacing between adjacent points in the Brillouin zone being $\delta k=0.1$
(left) and $\delta k=0.025$ (right).

###  III.6 Finite system size effects in disordered systems

Another manifestation of the difficulties associated with finite system sizes
is in calculating the topological invariant and band gaps for disordered
systems. In Fig. S6 we show the variation and (slow) convergence of both the
topological invariant, $N_{xy}$, as well as the bulk and edge band gaps as a
function of system size. Here, the system contains some disorder with a fixed
strength, $W/v=4.8$, with $W$ defined in the Main Text. The underlying system
possess $C_{4v}$ symmetry, with $W_{1}/v=1$ and $w_{2}/v=4$, while the added
disorder breaks all of the system’s symmetries (including $C_{4}$) except for
chiral symmetry. Only a single realization of the disorder is shown here at
each system size.

Figure S6: Dependence of the topological invariant $N_{xy}$, the bulk band
gap, and the two edge band gaps on the number of unit cells per side in a
square lattice. The underlying system has $C_{4v}$ symmetry, while the added
disorder breaks all of the symmetries (including $C_{4}$) except for chiral
symmetry.

##  IV Bulk and boundary energy gap closing at phase transitions

As shown in the Main Text, the phase transitions between phases of different
multipole chiral numbers close the gap at either the boundary or the bulk.
This has to be the case since changing the winding number across the phase
transition goes hand in hand with the delocalization of the topological corner
states, which necesitates extended gapless channels. The existence of these
gapless extended channels is necessary because the hybridization of
topological states away from zero energy can only occur by fusing pairs of
topological states having opposite chiral charge, and such pairs do not exist
within any single corner.

This is illustrated in Fig. S7. At the phase transition critical points, the
bulk (a) or edges (b,c) close the energy gap, allowing the topological corner
states to delocalize and hybridize away from zero energy. In (a), a $C_{4}$
symmetric lattice within a bulk-obstructed phase transition hybridizes corner
states via bulk low energy channels, reducing the bulk invariant by 1. In (b),
a $C_{4}$ symmetric boundary-obstructed phase transition hybridize the corner
states via the boundary, reducing the bulk invariant by 2. In (c), a $C_{2}$
symmetric boundary obstructed transition reduces the bulk invariant by 1.

Figure S7: Schematics that illustrates the hybridization of topological corner
states during phase transitions. Red and blue dots represent topological
corner states of opposite chiral charge. Topological corner states hybridize
only if they have opposite chiral charge and bulk (a) or edge (b,c) low energy
channels are available for their delocalization.

## References

  * Note [1] $G_{k}$ becomes unitary only in the thermodynamic limit, where the spectrum $k$ is continuous.
  * Lin  _et al._ [2021] L. Lin, Y. Ke, and C. Lee, Real-space representation of the winding number for a one-dimensional chiral-symmetric topological insulator, Phys. Rev. B 103, 224208 (2021).
  * Resta [1998] R. Resta, Quantum-mechanical position operator in extended systems, Phys. Rev. Lett. 80, 1800 (1998).
  * Teo and Kane [2010] J. C. Y. Teo and C. L. Kane, Topological defects and gapless modes in insulators and superconductors, Phys. Rev. B 82, 115120 (2010).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Topological quantum chemistry, Nature 547, 298 (2017).
  * Cano  _et al._ [2018] J. Cano, B. Bradlyn, Z. Wang, L. Elcoro, M. G. Vergniory, C. Felser, M. I. Aroyo, and B. A. Bernevig, Building blocks of topological quantum chemistry: Elementary band representations, Phys. Rev. B 97, 035139 (2018).


11institutetext:  Department of Physics, Chungbuk National University,
Cheongju 28644, Republic of Korea
11email<EMAIL_ADDRESS>22institutetext: Korea Astronomy and
Space Science Institute, Daejon 34055, Republic of Korea 33institutetext: Max
Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany
44institutetext: Department of Astronomy, The Ohio State University, 140 W.
18th Ave., Columbus, OH 43210, USA 55institutetext: University of Canterbury,
Department of Physics and Astronomy, Private Bag 4800, Christchurch 8020, New
Zealand 66institutetext: Center for Astrophysics $|$ Harvard & Smithsonian 60
Garden St., Cambridge, MA 02138, USA 77institutetext: Department of Particle
Physics and Astrophysics, Weizmann Institute of Science, Rehovot 76100, Israel
88institutetext: Department of Astronomy and Tsinghua Centre for Astrophysics,
Tsinghua University, Beijing 100084, China 99institutetext: School of Space
Research, Kyung Hee University, Yongin, Kyeonggi 17104, Republic of Korea
1010institutetext: Korea University of Science and Technology, 217 Gajeong-ro,
Yuseong-gu, Daejeon, 34113, Republic of Korea

# KMT-2021-BLG-1122L: The first microlensing triple stellar system

Cheongho Han 0101 Youn Kil Jung 0202 Andrew Gould 03030404 Doeon Kim 0101
Chung-Uk Lee
(Leading authors)
and
Michael D. Albrow 02020505 Sun-Ju Chung 0202 Kyu-Ha Hwang 0202 Hyoun-Woo Kim
0202 Yoon-Hyun Ryu 0202 In-Gu Shin 0606 Yossi Shvartzvald 0707 Hongjing Yang
0808 Jennifer C. Yee 0606 Weicheng Zang 06060808 Sang-Mok Cha 02020909 Dong-
Jin Kim 0202 Seung-Lee Kim 02021010 Dong-Joo Lee 0202 Yongseok Lee 02020909
Byeong-Gon Park 02021010 Richard W. Pogge
(The KMTNet Collaboration)
0404

(Received ; accepted)

###### Abstract

Aims.  We systematically inspect the microlensing data acquired by the KMTNet
survey during the previous seasons in order to find anomalous lensing events
for which the anomalies in the lensing light curves cannot be explained by the
usual binary-lens or binary-source interpretations.

Methods.  From the inspection, we find that interpreting the three lensing
events OGLE-2018-BLG-0584, KMT-2018-BLG-2119, and KMT-2021-BLG-1122 requires
four-body (lens+source) models, in which either both the lens and source are
binaries (2L2S event) or the lens is a triple system (3L1S event). Following
the analyses of the 2L2S events presented in Han et al. (2023), here we
present the 3L1S analysis of the KMT-2021-BLG-1122.

Results.  It is found that the lens of the event KMT-2021-BLG-1122 is composed
of three masses, in which the projected separations (normalized to the angular
Einstein radius) and mass ratios between the lens companions and the primary
are $(s_{2},q_{2})\sim(1.4,0.53)$ and $(s_{3},q_{3})\sim(1.6,0.24)$. By
conducting a Bayesian analysis, we estimate that the masses of the individual
lens components are
$(M_{1},M_{2},M_{3})\sim(0.47\,M_{\odot},0.24\,M_{\odot},0.11\,M_{\odot})$.
The companions are separated in projection from the primary by
$(a_{\perp,2},a_{\perp,3})\sim(3.5,4.0)$ AU. The lens of KMT-2018-BLG-2119 is
the first triple stellar system detected via microlensing.

###### Key Words.:

Gravitational lensing: micro – (Stars:) binaries: general

##  1 Introduction

Light curves of gravitational microlensing events often exhibit deviations
from the smooth and symmetric form of the event involved with a single lens
mass and a single source star (Paczyński, 1986). The most common cause of such
deviations is the binarity of the lens, 2L1S event (Mao & Paczyński, 1991), or
the source, 1L2S event (Griest & Hu, 1993; Han & Gould, 1997). Anomalous
lensing events comprise about 10% of more than 3000 lensing events that are
being annually detected from the three currently working lensing surveys
conducted by the Optical Gravitational Lensing Experiment (OGLE: Udalski et
al., 2015), Microlensing Observations in Astrophysics (MOA: Bond et al.,
2001), and Korea Microlensing Telescope Network (KMTNet: Kim et al., 2016)
groups. The prime goal of these survey experiments is finding extrasolar
planets, which are detected via anomalous signals in lensing light curves.
Planetary microlensing signals appear in various shapes and often confused
with anomalies produced by other causes. In order to sort out planetary
signals from those of other origins, therefore, it is essential to analyze all
anomalous events.

In the current microlensing experiments, anomalous events are being analyzed
almost in real time with the progress of events. This became possible with the
implementation of early warning systems: the OGLE early warning system
(Udalski et al., 1994), MOA real-time analysis system (Bond et al., 2001), and
KMTNet EventFinder AlertFinder (Kim et al., 2018b) system.111In addition, the
KMTNet EventFinder system (Kim et al., 2018a) carries out end-of-season
analyses that recover many events that were missed for one reason or another
by AlertFinder (Kim et al., 2018b) system. These systems enable one to detect
lensing events in their early phases and to promptly identify anomalies
appearing in lensing light curves. Analyses of the events with identified
anomalies are conducted by multiple modelers of the individual survey groups
as the anomalies proceed, and models of the events found from these analyses
are circulated or posted on web pages222For example, the model page of lensing
events maintained by Cheongho Han
(http://astroph.chungbuk.ac.kr/$\sim$cheongho). to inform researchers in the
microlensing community of the nature of the anomalies. In most cases,
anomalous events are interpreted with a binary-lens (2L1S) or a binary-source
(1L2S) model.

For a minor fraction of anomalous events, it is known that the anomalies
cannot be explained by the usual 2L1S or 1L2S interpretation. Han et al.
(2023), hereafter Paper I, conducted a systematic investigation of events
found from the KMTNet survey during the previous seasons in search for
anomalous events for which no plausible models had been presented. From this
investigation, they found that the anomalies appearing in the light curves of
the events OGLE-2018-BLG-0584 and KMT-2018-BLG-2119 required a four-body
(lens+source) model, in which both the lens and source are binaries (2L2S
event). In this paper, we present the analysis of another four-body lensing
event KMT-2021-BLG-1122. Unlike the two 2L2S events presented in Paper I, the
lens system of KMT-2021-BLG-1122 is composed of three masses (3L1S).

We present the analysis of the event KMT-2021-BLG-1122 according to the
following organization. In Sect. 2, we depict the observations conducted to
obtain the photometric data of the event, and describe the characteristic
features of the anomaly appearing in the light curve of the event. In Sect. 3,
we begin by describing the lensing parameters used in the modeling conducted
under various configurations of the lens and source system, and we then
describe detail of the analyses conducted under the individual configurations
in the following subsections: 2L1S model in Sect. 3.1, 2L2S model in Sect.
3.2, and 3L1S model in Sect. 3.3. In Sect. 4, we specify the source star of
the event and estimate the angular Einstein radius of the lens system. In
Sect. 5, we estimate the physical parameters including the masses of the
individual lens components, the distance to the lens system, and the projected
separations among the lens components. We summarize the results from the
analysis and conclude in Sect. 6.

##  2 Lensing light curve and observations

The lensing event KMT-2021-BLG-1122 was detected solely by the KMTNet survey.
The event occurred on a faint source star, with a baseline magnitude of
$I_{\rm base}=20.36$, located toward the Galactic bulge field with equatorial
coordinates $({\rm RA},{\rm DEC})=$ (17:35:51.10, $-$28:26:47.22), which
correspond to the Galactic coordinates
$(l,b)=(-0^{\circ}\hskip-2.0pt.720,2^{\circ}\hskip-2.0pt.076)$. The source
lies in the KMTNet BLG14 field, toward which observations were conducted with
a 1 hr cadence. The event was identified by the KMTNet AlertFinder system and
an alert was issued on 2021 June 3, which corresponds to the abridged
Heliocentric Julian Date of ${\rm HJD}^{\prime}\equiv{\rm HJD}-2450000\sim
9368.5$.

Observations of the event were carried out with the use of the three 1.6 m
telescopes operated by the KMTNet group. The telescopes are distributed in the
three continents of the Southern Hemisphere for continuous coverage of lensing
events, and the sites of the individual telescopes are the Siding Spring
Observatory in Australia (KMTA), the Cerro Tololo Interamerican Observatory in
Chile (KMTC), and the South African Astronomical Observatory in South Africa
(KMTS). The camera mounted on each telescope is composed of four $9{\rm
k}\times 9{\rm k}$ CCD chips, which in combination yield a 4 deg2 field of
view.

Images containing the source were mainly taken in the $I$ band, and about one
tenth of images were obtained in the $V$ band. Reductions of the images and
photometry of the event were done using the KMTNet pipeline constructed based
on the pySIS code of Albrow et al. (2009), and additional photometry was
carried out for a subset of the KMTC data set using the pyDIA code of Albrow
(2017). The analysis of the event was done based on the $I$-band light curve
constructed from the pySIS reduction, and the pyDIA photometry data were used
for the source color measurement. We discuss the detailed procedure of the
source color measurement in Sect. 4. According to recipe described in Yee et
al. (2012), we readjusted the error bars of the data determined by the
automated pipeline in order that the errors were consistent with the scatter
of data and $\chi^{2}$ per degree of freedom (dof) for each data set became
unity.

Figure 1:  Light curve of KMT-2021-BLG-1122 constructed with the combined data
obtained from the three KMTNet telescopes (KMTA, KMTC, and KMTS). The curve
drawn over the data points is the best-fit model found under the 2L1S
interpretation. The lower panel displays the whole view, and the upper panel
shows the enlarged view of the anomaly region. The epochs marked by $t_{1}$,
$t_{2}$, and $t_{3}$ indicate the three major anomaly features.  Figure 2:
Enlarged view around the central anomaly region of the light curve. Three
model curves obtained under the 3L1S, 2L2S, and 2L1S lens-system
configurations are drawn over the data points, and the three lower panels show
the residuals from the individual models. The two insets in the top panel show
the lens-system configurations of the 2L1S and 2L2S models. In each inset, the
red figure is the caustic and the line with an arrow represents the source
trajectory. For the 2L2S configuration, there are two source trajectories, in
which the black and blue ones represent the trajectories of the primary and
secondary source stars, respectively.

Figure 1 shows the lensing light curve of KMT-2021-BLG-1122 constructed with
the combined $I$-band pySIS photometry data obtained from the three KMTNet
telescopes. It shows a strong anomaly relative to a 1L1S form, and the anomaly
displays a complex pattern that is characterized by the three strong features
appearing at ${\rm HJD}^{\prime}\sim 9369.7$ ($t_{1}$), $\sim 9371.7$
($t_{2}$), and $\sim 9373.3$ ($t_{3}$). From the characteristic profiles
described in Schneider & Weiss (1986), the anomalies at around $t_{1}$ and
$t_{3}$ appear to be a pair of caustic-crossing features produced when a
source enters and exits a caustic, respectively. On the other hand, the
profile at around $t_{2}$, which is approximately symmetric with respect to
the peak, appears to be produced when a source approaches close to or crosses
over the cusp of a caustic (Schneider & Weiss, 1987).

In general, the light curve profile of a 2L1S event in the region between a
pair caustic-crossing spikes exhibits a ”U”-shape pattern. Occasionally,
deviations from this U-shape pattern can arise if a source sweeps one fold of
a caustic, for example, OGLE-2016-BLG-0890 (Han et al., 2022c), but, in this
case, the resulting deviation is smooth and thus differs from the sharp
pattern of the observed anomaly at around $t_{2}$. This suggests that a
different explanation is needed to explain the origin of the anomaly.

##  3 Interpreting the anomaly

In order to explain the anomalies in the lensing light curve, we conducted a
series of modeling under various interpretations of the lens-system
configuration. We tested three models under the 2L1S, 2L2S, and 3L1S
configurations. We exclude 1L2S and 1L3S models because the observed light
curve exhibits obvious caustic-crossing features, that is, those around
$t_{1}$ and $t_{3}$, while 1L2S and 1L3S events do not involve
caustics.333There has been only one 1L3S microlensing event published to date:
OGLE-2015-BLG-1459 (Hwang et al., 2018).

In the modeling under each tested interpretation, we search for a lensing
solution, representing a set of lensing parameters that characterize the
observed lensing light curve. The basic lensing parameters of a 1L1S event
include $(t_{0},u_{0},t_{\rm E})$, which denote the time of the closest source
approach to the lens, the source-lens separation (impact parameter) at
$t_{0}$, and the Einstein time scale, respectively. The Einstein time scale
represents the time required for the source to transit the angular Einstein
radius of the lens, $\theta_{\rm E}$, and the length of $u_{0}$ is scaled to
$\theta_{\rm E}$.

A 2L1S system corresponds to the case in which the lens contains an extra
component compared to the 1L1S system. Consideration of the extra lens
component requires one to include additional parameters in modeling. These
additional parameters are $(s,q,\alpha,\rho)$, and they denote the projected
separation (normalized to $\theta_{\rm E}$), and mass ratio between the lens
components ($M_{1}$ and $M_{2}$), the angle (source trajectory angle) between
the relative lens-source proper motion, $\mu$, and the axis connecting $M_{1}$
and $M_{2}$, and the ratio between the angular radius of the source,
$\theta_{*}$, to the Einstein radius, that is, $\rho=\theta_{*}/\theta_{\rm
E}$ (normalized source radius), respectively. The normalized source radius is
included in modeling because a 2L1S event usually involves caustic crossings,
during which lensing magnifications are affected by finite-source effects
(Bennett & Rhie, 1996). For the computation of finite-source magnifications,
we utilize the map-making method of Dong et al. (2008).

For the modeling of a 2L2S system, it is required to include extra parameters
in addition to those of the 2L1S model in order to describe an extra source.
These extra parameters include $(t_{0,2},u_{0,2},\rho_{2},q_{F})$, where the
first three denote the closest approach time, impact parameter, and normalized
radius of the source companion, and the last parameter indicates the flux
ratio between the companion ($S_{2}$) and primary ($S_{1}$) of the source. In
the 2L2S modeling, we use the notations $(t_{0,1},u_{0,1},\rho_{1})$ to denote
the parameters related to $S_{1}$.

A 3L1S system has an extra lens component ($M_{3}$) compared to the 2L1S
system, and this requires one to include additional parameters in modeling.
These parameters include $(s_{3},q_{3},\psi)$, which denote the separation and
mass ratio between $M_{1}$ and $M_{3}$, and the orientation angle of $M_{3}$
as measured from the $M_{1}$–$M_{2}$ axis with a center at the position of
$M_{1}$, respectively. In order to distinguish the parameters describing
$M_{2}$ from those describing $M_{3}$, we use the notations $(s_{2},q_{2})$ to
denote the separation and mass ratio of $M_{2}$. The lensing parameters for
the 1L1S, 2L1S, 2L2S, and 3L1S models are summarized in Table 1 of Paper I.

###  3.1 2L1S interpretation

We started with modeling the light curve under the 2L1S interpretation of the
anomalies. In this modeling, we searched for the binary parameters $s$ and $q$
via a grid approach, while we seek the other parameters via a downhill
approach using a Markov Chain Monte Carlo (MCMC) algorithm. We first
constructed a $\Delta\chi^{2}$ map on the plane of the grid parameters,
identified local minima, and then refined the individual minima by allowing
all lensing parameters to vary. From the 2L1S modeling, we found no lensing
solution that could describe all the anomaly features, and this explains why
the event had been left without a plausible model.

Although the light curve could not be explained by a 2L1S model, we found a
solution that could approximately describe the caustic-crossing features at
around $t_{1}$ and $t_{3}$. The model curve of this solution, found from
modeling the light curve with the exclusion of the data around $t_{2}$, is
drawn over the data points in Figures 1 and 2, and the corresponding
configuration of the lens system is presented in the upper inset of the top
panel of Figure 2. The binary lensing parameters of the solution are
$(s,q)\sim(1.69,0.43)$, which results in a single resonant caustic elongated
along the binary axis. According to this solution, the source entered the
caustic by crossing the caustic fold lying just above the left side on-axis
cusp, passed through the inner region of the caustic, and then exited the
caustic by crossing the upper fold of the caustic. The caustic entrance and
exit resulted in the anomaly features matching well the observed ones around
$t_{1}$ and $t_{3}$, respectively.

Table 1: Lensing parameters of the 2L2S solutions

Parameter | Static | Orbital
---|---|---
$\chi^{2}$/dof | 1729.1/1426 | 1453.6/1429
$t_{0,1}$ (HJD′) | $9371.208\pm 0.026$ | $9370.681\pm 0.007$
$u_{0,2}$ | $0.057\pm 0.002$ | $0.044\pm 0.001$
$t_{0,2}$ (HJD′) | $9370.810\pm 0.024$ | –
$u_{0,2}$ | $0.194\pm 0.004$ | –
$t_{\rm E}$ (days) | $14.36\pm 0.31$ | $14.87\pm 0.11$
$s$ | $1.723\pm 0.020$ | $1.447\pm 0.006$
$q$ | $0.649\pm 0.043$ | $0.261\pm 0.002$
$\alpha$ (rad) | $2.615\pm 0.009$ | $2.415\pm 0.005$
$s_{s}$ | – | $0.079\pm 0.001$
$q_{s}$ | – | $0.103\pm 0.001$
$\phi$ (rad) | – | $1.326\pm 0.012$
$ds_{s}/dt$ (yr-1) | – | $-0.064\pm 0.031$
$d\phi/dt$ (deg/day) | – | $53.79\pm 0.30$
$\rho_{1}$ (10-3) | $3.29\pm 0.17$ | $4.58\pm 0.05$
$\rho_{2}$ (10-3) | $2.23\pm 0.79$ | $0.10\pm 0.01$
$q_{F}$ | $0.71\pm 0.05$ | $0.15\pm 0.01$

444 ${\rm HJD}^{\prime}={\rm HJD}-2450000$.

###  3.2 2L2S interpretation

The fact that the two caustic-crossing features can be described by a 2L1S
model suggests that the other anomaly feature around $t_{2}$ may be explained
with the introduction of an extra source (2L2S) or an extra lens component
(3L1S) to the 2L1S configuration. In this section, we present the analysis
conducted under the 2L2S interpretation.

The 2L2S modeling was carried out based on the 2L1S solution. With the initial
parameters obtained from the 2L1S modeling, we searched for a 2L2S solution by
testing various trajectories of the second source considering the location and
magnitude of the anomaly around $t_{2}$. We first check the case in which the
position of the source companion with respect to the primary source does not
vary: ”static 2L2S model”. The model curve of the best-fit static 2L2S
solution is presented in Figure 2, and the lensing parameters are listed in
Table 1. The lens-system configuration of the solution is shown in the middle
inset of the top panel. The caustic configuration is similar to that of the
2L1S solution except that there are two source trajectories of $S_{1}$ and
$S_{2}$, which are marked by black and blue arrowed lines, respectively.
According to the solution, the source companion approached the lens slightly
earlier than the primary source with an impact parameter greater than that of
the primary source. The second source passed the tip of the upper left caustic
cusp, and this produced the anomaly feature around $t_{2}$. Although the
static 2L2S model can give rise to a sharp anomaly features appearing in the
inner region between the caustic spikes, it was found that the residual from
the 2L2S model was substantial, as shown in the second lower panel of Figure
2.

Figure 3:  Map of $\Delta\chi^{2}$ on the $\log\leavevmode\nobreak\
s_{3}$–$\log\leavevmode\nobreak\ q_{3}$ parameter plane obtained from the
initial grid searches for $(s_{3},q_{3})$ parameters. Color coding is set to
represent points with $\leq 1n\sigma$ (red), $\leq 2n\sigma$ (yellow), $\leq
3n\sigma$ (green), $\leq 4n\sigma$ (cyan), $\leq 5n\sigma$ (blue), and $\leq
6n\sigma$ (purple), where $n=9$. The upper panel shows the map around the
minima constructed with denser grids.

We further check whether the anomalies can be explained by considering the
orbital motion of the source: ”orbital 2L2S model”. We consider the source
orbital motion by introducing five extra parameters of ($s_{s}$, $q_{s}$,
$\phi$, $ds_{s}/dt$, $d\phi/dt$), where $s_{s}$ and $q_{s}$ represent the
normalized separations and mass ratio between the source components,
respectively, $\phi$ is the orientation angle of $S_{2}$ with respect to
$S_{1}$, and $(ds_{s}/dt,d\phi/dt)$ represent the change rates of the $s_{s}$
and $\phi$ induced by the source orbital motion, respectively. It is found
that the consideration of the source orbital motion substantially improves the
fit, by $\Delta\chi^{2}=275.5$, with respect to the static model. The model
curve of the best-fit orbital 2L2S solution is drawn in the top panel of
Figure 2, residual is shown in the panel labeled ”2L2S (orbit)”, and the
lensing parameters are listed in Table 1. According to the model, the anomaly
features around $t_{1}$ and $t_{3}$ were produced by the primary source as is
in the static case, and the feature around $t_{2}$ was produced by the
secondary source. We note that the secondary source crossed the caustic two
more times at ${\rm HJD}^{\prime}\sim 9368.7$ ($t_{0}$) and $\sim 9372.7$, and
the former crossing is not covered by the data and the latter crossing feature
is difficult to identify in the data. We note that the source radius ratio
$\rho_{2}/\rho_{1}=1/46$ and the source flux ratio $q_{F}=0.15$ would lead to
inconsistent estimates of $\theta_{\rm E}$, which would argue against the
physical plausibility of this solution. However, we do not pursue this issue
because this solution will prove to be heavily disfavored by $\chi^{2}$. See
Sect. 3.3.

Figure 4:  Lens-system configuration of the 3L1S model. The main panel shows
the central anomaly region, and the inset displays the whole view of the lens
system, where blue dots marked by $M_{1}$, $M_{2}$, and $M_{3}$ denote the
positions of the individual lens components. Grey curves around the caustic
represent the equi-magnification contours. The dotted circle centered at the
origin of the coordinates represents the Einstein ring of the lens system.
Table 2: Lensing parameters of the best-fit 3L1S solution

Parameter | Value
---|---
$\chi^{2}$/dof | 1399.8/1425
$t_{0}$ (HJD′) | $9370.226\pm 0.061$
$u_{0}$ | $0.126\pm 0.008$
$t_{\rm E}$ (days) | $14.74\pm 0.93$
$s_{2}$ | $1.386\pm 0.027$
$q_{2}$ | $0.526\pm 0.039$
$\alpha$ (rad) | $2.292\pm 0.013$
$s_{3}$ | $1.601\pm 0.043$
$q_{3}$ | $0.241\pm 0.027$
$\psi$ (rad) | $1.383\pm 0.010$
$\rho$ ($10^{-3}$) | $2.50\pm 0.21$

555 ${\rm HJD}^{\prime}={\rm HJD}-2450000$.

###  3.3 3L1S interpretation

We additionally searched for a model under the 3L1S interpretation. Similar to
the case of the 2L2S modeling, the 3L1S modeling was done based on the 2L1S
solution, from which we adopted the initial parameters of $(t_{0},u_{0},t_{\rm
E},s_{2},q_{2},\alpha,\rho)$. We then conducted preliminary searches for the
parameters related to $M_{3}$, that is, $(s_{3},q_{3},\psi)$, via a grid
approach, and finally refined the solution by letting all parameters vary.
Figure 3 shows the $\Delta\chi^{2}$ map on the $\log\leavevmode\nobreak\
s_{3}$–$\log\leavevmode\nobreak\ q_{3}$ parameter plane constructed from the
preliminary grid search. The map shows that there exists a unique $\chi^{2}$
minimum at $(\log\leavevmode\nobreak\ s_{3},\log\leavevmode\nobreak\
q_{3})\sim(0.2,-0.6)$.

It was found that the model identified under the 3L1S interpretation well
described all the observed anomaly features in the lensing light curve. The
lensing parameters of the best-fit 3L1S solution are listed in Table 2
together with the value of $\chi^{2}/{\rm dof}$. In Figure 2, we draw the
model curve (solid curve drawn over the data points) and residual (presented
in the panel labeled ”3L1S”) from the model. The parameters related to $M_{2}$
and $M_{3}$ are $(s_{2},q_{2})\sim(1.4,0.53)$ and
$(s_{3},q_{3})\sim(1.60,0.24)$, respectively, indicating that the primary of
the lens is accompanied by two lower-mass companions lying in the vicinity of
the Einstein ring. If the separations among the lens components were intrinsic
(3-dimensional), the three-body lens system would be dynamically unstable, and
this indicates that one or both companions are likely to be aligned by chance
with the primary. We found that it was difficult to detect microlens-parallax
effects (Gould, 1992) due to the relatively short time scale, $t_{\rm E}\sim
15$ days, of the event compared to the orbital period of Earth.

The lens-system configuration of the 3L1S model is shown in Figure 4, in which
the main panel shows the central anomaly region and the inset shows the whole
view of the lens system. The coordinates of the configuration is centered at
the effective position of the $M_{1}$–$M_{2}$ pair defined by Di Stefano & Mao
(1996) and An & Han (2002), and the grey curves encompassing the caustic
represent the equi-magnification contours. One finds that the caustic induced
by $M_{2}$ is similar to that of the 2L1S solution, but the 3L1S caustic and
the resulting magnification pattern around the caustic are distorted by
$M_{3}$, which makes the central part of the caustic nested and self-
intersected (Gaudi et al., 1998; Danĕk & Heyrovský, 2015a, b, 2019). According
to the 3L1S solution, the anomaly appearing around $t_{2}$ was produced by the
source approach close to the tip of the swallow-tail caustic part induced by
$M_{3}$.

Figure 5:  Cumulative distribution of the $\chi^{2}$ difference,
$\Delta\chi^{2}=\chi^{2}_{\rm 2L2S}-\chi^{2}_{\rm 3L1S}$, between the 2L2S
(orbital) and 3L1S models in the region of the major anomaly. The light curve
in the upper panel is shown to indicate the epochs of major $\chi^{2}$
differences.

From the comparison of the 3L1S and 2L2S (orbital) models, we find that the
anomaly features in the lensing light curve are better explained by the 3L1S
interpretation than by the 2L2S interpretation. In Figure 5, we present the
cumulative distribution of the $\chi^{2}$ difference,
$\Delta\chi^{2}=\chi^{2}_{\rm 2L2S}-\chi^{2}_{\rm 3L1S}$, between the two
models in the region of the anomaly for the detailed comparison of the two
models. The major differences occur around $t_{0}$ and $t_{1}$, at around
which $S_{2}$ and $S_{1}$ first crossed the caustic according to the 2L2S
model, respectively. Besides these regions, the 3L1S model well explains the
KMTA data point at the epoch of ${\rm HJD}^{\prime}=9373.272$, while the
residual of this data point from the 2L2S model, $\Delta I=0.656$ mag, is very
big. As a whole, the 3L1S model yields a better fit than the 2L2S model by
$\Delta\chi^{2}=53.8$, despite the fact that the 3L1S model has 4 fewer
parameters.

##  4 Source star and Einstein radius

We specify the source of the event not only for the full characterization of
the event but also for the estimation of the angular Einstein radius. The
Einstein radius is determined from the angular radius of the source,
$\theta_{*}$, by the relation

$\theta_{\rm E}={\theta_{*}\over\rho},$ (1)

where $\theta_{*}$ is deduced from the source type, and the normalized source
radius $\rho$ is measured from the modeling by analyzing the caustic-crossing
parts of the light curve.

We determined the source type by measuring the extinction- and reddening-
corrected color and brightness, $(V-I,I)_{0,S}$. For this, we utilize the Yoo
et al. (2004) method, in which the source is located in the instrumental
color-magnitude diagram (CMD), and then its color and magnitude are calibrated
using the centroid of the red giant clump (RGC) as a reference. The RGC
centroid can be used as a reference for calibration because its de-reddened
color and magnitude are known.

Figure 6 shows the position of the source in the instrumental CMD constructed
from the pyDIA photometry of stars lying around the source. We note that the
CMD shows only bright stars due to the severe extinction, $A_{I}=3.3$, toward
the field. The measured instrumental color and magnitude of the source are

$(V-I,I)_{S}=(3.501\pm 0.056,21.080\pm 0.037),$ (2)

where the $I$ and $V$-band magnitudes were estimated by regressing the KMTC
light curve data of the individual passbands processed using the same pyDIA
code with respect to the model. By measuring the offsets in the color and
magnitude, $\Delta(V-I,I)$, of the source from the RGC centroid, with
$(V-I,I)_{\rm RGC}=(3.669,18.025)$, together with the known de-reddened color
and magnitude of the RGC centroid, $(V-I,I)_{0,{\rm RGC}}=(1.060,14.377)$
(Bensby et al., 2013; Nataf et al., 2013), the de-reddened source color and
magnitude were measured as

$\displaystyle{(V-I,I)_{0,S}}$ | $\displaystyle{{}=(V-I,I)_{0,{\rm RGC}}+\Delta(V-I,I)}$
---|---
| $\displaystyle{{}=(0.892\pm 0.056,17.432\pm 0.037).}$
(3)

The estimated color and magnitude indicate that the source is a subgiant star
of an early K spectral type.

Figure 6:  Position of the source (blue dot) with respect to the centroid of
red giant clump (RGC, red dot) in the instrumental color-magnitude diagram of
stars lying in the vicinity of the source.

The angular source radius was deduced from the measured source color and
magnitude. For this, we first converted $V-I$ color into $V-K$ color using the
color-color relation of Bessell & Brett (1988), and then estimated the angular
source radius using the relation between $(V-K,V)$ and $\theta_{*}$ of
Kervella et al. (2004). This procedure yielded the angular radii of the source
and Einstein ring of

$\theta_{*}=(1.26\pm 0.11)\leavevmode\nobreak\ \mu{\rm as}.$ (4)

and

$\theta_{\rm E}={\theta_{*}\over\rho}=(0.50\pm 0.06)\leavevmode\nobreak\ {\rm
mas},$ (5)

respectively. Together with the measured event time scale, the relative lens-
source proper motion was estimated as

$\mu={\theta_{\rm E}\over t_{\rm E}}=(12.49\pm 1.53)\leavevmode\nobreak\ {\rm
mas}\leavevmode\nobreak\ {\rm yr}^{-1}.$ (6)

##  5 Physical lens parameters

In this section, we estimate the physical lens parameters of the mass $M$ and
distance $D_{\rm L}$ to the lens. These parameters can be sometimes
constrained by measuring lensing observables of $(t_{\rm E},\theta_{\rm
E},\pi_{\rm E})$, where $\pi_{\rm E}$ denotes the microlens parallax. The
first two observables are related to the physical parameters by

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

where $\kappa=4G/(c^{2}{\rm AU})\simeq 8.14\leavevmode\nobreak\ {\rm
mas}/M_{\odot}$ and $\pi_{\rm rel}=\pi_{\rm L}-\pi_{\rm S}={\rm AU}(1/D_{\rm
L}-1/D_{\rm S})$ represents the relative lens-source parallax. With the
additional measurement of the observable $\pi_{\rm E}$, the physical
parameters can be uniquely constrained by the relations (Gould, 2000)

$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}}.$ (8)

For KMT-2021-BLG-1122, the observables $t_{\rm E}$ and $\theta_{\rm E}$ were
securely measured, but the other observable $\pi_{\rm E}$ could not be
measured, and thus we estimate the physical lens parameters by conducting a
Bayesian analysis using a Galactic model together with the constraints
provided by the measured observables.

In the first step of the Bayesian analysis, a large number of artificial
lensing events were produced from a Monte Carlo simulation. For each simulated
event, we derived the locations of the lens and source and their transverse
velocity from a Galactic model, and assigned the lens mass from a mass-
function model. In the Monte Carlo simulation, we adopted the Galactic model
of Jung et al. (2021), whose detailed description are given therein. With the
simulated events, we then constructed the posteriors of the lens mass and
distance by giving a weight to each simulated event of

$w_{i}=\exp(-\chi^{2}_{i}/2);\qquad\chi_{i}^{2}={(t_{{\rm E},i}-t_{\rm
E})^{2}\over[\sigma(t_{\rm E})]^{2}}+{(\theta_{{\rm E},i}-\theta_{\rm
E})^{2}\over[\sigma(\theta_{\rm E})]^{2}},$ (9)

where $(t_{\rm E},\theta_{\rm E})$ and $[\sigma(t_{\rm E}),\sigma(\theta_{\rm
E})]$ represent the measured values of the observables and their
uncertainties, respectively, and $(t_{{\rm E},i},\theta_{{\rm E},i})$ indicate
the observables of each simulated event.

Figure 7:  Posteriors of the primary lens mass (top panel) and distance to the
lens (middle panel) and source (bottom panel). In each panel, the blue and red
curves represent the contributions by the disk and bulge lens populations,
respectively, and the black curve is sum of the two lens populations. The
solid vertical line indicates the median value, and the two dotted lines
represent the $1\sigma$ range of the distribution.

Figure 7 shows the posteriors of the primary lens mass (top panel), and
distances to the lens (middle panel) and source (bottom panel). We
additionally present the $D_{\rm S}$ posterior to show the relative location
of the lens and source. For each posterior, we mark the contributions by the
disk and bulge lens populations with curves drawn in blue and red,
respectively. The masses of the individual lens components estimated from the
Bayesian analysis are

$M_{1}=0.47^{+0.24}_{-0.33}\leavevmode\nobreak\ M_{\odot},$ (10)
$M_{2}=0.24^{+0.13}_{-0.17}\leavevmode\nobreak\ M_{\odot},$ (11)

and

$M_{3}=0.11^{+0.06}_{-0.08}\leavevmode\nobreak\ M_{\odot}.$ (12)

The estimated masses of all the lens components are in the mass regime of M
dwarfs, and thus the lens is a triple system composed of three low-mass stars.
The estimated distance to the lens system is

$D_{\rm L}=6.57^{+1.34}_{-0.94}\leavevmode\nobreak\ {\rm kpc},$ (13)

and the probabilities for the lens to be in the disk and bulge are 33% and
67%, respecively. For each parameter, we present the median as a
representative value and the lower and upper limits are estimated as the 16%
and 84% ranges of the posterior distribution. The projected separations of the
secondary and tertiary lens components from the primary are

$a_{\perp,2}=s_{2}D_{\rm L}\theta_{\rm
E}=3.47^{+0.71}_{-0.50}\leavevmode\nobreak\ {\rm AU},$ (14)

and

$a_{\perp,3}=s_{3}D_{\rm L}\theta_{\rm
E}=4.01^{+0.82}_{-0.57}\leavevmode\nobreak\ {\rm AU},$ (15)

respectively. As noted, the fact that the projected separations $a_{\perp,2}$
and $a_{\perp,3}$ are similar to each other suggests that one or both lens
companions are aligned by chance with the primary.

We note that KMT-2021-BLG-1122L is the first triple stellar system detected
via microlensing. In total, there exist 12 confirmed 3L1S events:
OGLE-2006-BLG-109 (Gaudi et al., 2008; Bennett et al., 2010),
OGLE-2012-BLG-0026 (Han et al., 2013), OGLE-2018-BLG-1011 (Han et al., 2019),
OGLE-2019-BLG-0468 (Han et al., 2022d), OGLE-2021-BLG-1077 (Han et al.,
2022a), OGLE-2006-BLG-284 (Bennett et al., 2020), OGLE-2007-BLG-349 (Bennett
et al., 2016), OGLE-2008-BLG-092 (Poleski et al., 2014), OGLE-2013-BLG-0341
(Gould et al., 2014), OGLE-2018-BLG-1700 (Han et al., 2020b),
KMT-2019-BLG-1715 (Han et al., 2021a), and KMT-2020-BLG-0414 (Zang et al.,
2021).666 Besides these confirmed cases, there are 6 candidate triple-lensing
events, including OGLE-2014-BLG-1722 (Suzuki et al., 2018), OGLE-2018-BLG-0532
(Ryu et al., 2020), KMT-2019-BLG-1953 (Han et al., 2022a), KMT-2019-BLG-0304
(Han et al., 2021b), OGLE-2019-BLG-1470 (Kuang et al., 2022), and
KMT-2021-BLG-0240 (Han et al., 2022b). For these events, 3L1S interpretations
yielded best-fit models, but the signals of the tertiary lens components were
not firmly confirmed either due to the subtlety of the signals or the
degeneracy with other interpretations. Under the 3L1S interpretations of these
events, the lenses belong to either the category of two-planet systems or the
category of planets in binaries similar to the cases of the confirmed 3L1S
events. Among them, the lenses of the first 5 events are two-planet systems,
and the lenses of the other 7 events are binary-star systems possessing
planets. Therefore, the common trait of the previously discovered microlensing
triple systems is that at least one of the lens components is a planet. Two
planets with similar separations from a host can be dynamically stable if the
planets are in mean-motion resonance (Madsen & Zhu, 2019). A planet in a
binary system can be dynamically stable if the planet orbits one star of a
wide binary stellar system or the barycenter of a closely-spaced binary
system. Microlensing detection of a triple stellar system is difficult because
a system with similar masses and intrinsic separations among the components
would be dynamically unstable. One way for such a system to be detected via
microlensing is that the companions of the triple system are closely aligned
with the primary so that they lie around the Einstein ring of the primary.
This alignment requires projection effect of a low probability, and this
explains the relative rareness of triple stellar systems detected by
microlensing.

##  6 Summary and conclusion

We presented the analysis of the microlensing event KMT-2021-BLG-1122, which
was investigated as a part of the project to analyze anomalous lensing events
in the previous KMTNet data with no suggested plausible models. We confirmed
that the light curve, characterized by three major anomaly features, of the
event could not be explained with the usual 2L1S or 1L2S interpretations, but
it was found that a 2L1S solution obtained from the modeling with the
exclusion of the data around the second anomaly feature could approximately
describe the other two caustic-crossing features.

We found that all the anomaly features could be explained by introducing a
tertiary lens component, and thus the lens is a triple system. The masses of
the individual lens components estimated from the Bayesian analysis are
$(M_{1},M_{2},M_{3})\sim(0.47\leavevmode\nobreak\
M_{\odot},0.24\leavevmode\nobreak\ M_{\odot},0.11\leavevmode\nobreak\
M_{\odot})$, and the companions are separated in projection from the primary
by $(a_{\perp,2},a_{\perp,3})\sim(3.5,4.0)$ AU. The estimated masses of all
the lens components are in the mass regime of M dwarfs, and thus the lens is a
triple system composed of three low-mass stars. This is the first triple
stellar-lens system detected via microlensing.

###### Acknowledgements.

Work by C.H. was supported by the grants of National Research Foundation of
Korea (2020R1A4A2002885 and 2019R1A2C2085965). J.C.Y. acknowledges the support
from NSF Grant No. AST-2108414. W.Z. and H.Y. acknowledge the support by the
National Science Foundation of China (Grant No. 12133005). Y.S. acknowledges
support from BSF Grant No. 2020740. 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

  * Albrow (2017) Albrow, M. 2017, MichaelDAlbrow/pyDIA: Initial Release on Github,Versionv1.0.0, Zenodo, doi:10.5281/zenodo.268049
  * Albrow et al. (2009) Albrow, M., Horne, K., Bramich, D. M., et al. 2009, MNRAS, 397, 2099
  * An & Han (2002) An, J. H., & Han, C. 2002, ApJ, 573, 351
  * Beaulieu et al. (2016) Beaulieu, J.-P., Bennett, D. P., Batista, V., et al. 2016, ApJ, 824, 83
  * Bennett & Rhie (1996) Bennett, D. P., & Rhie, S. H. 1996, ApJ, 472, 660
  * Bennett et al. (2010) Bennett, D. P., Rhie, S. H., Nikolaev, S., et al. 2010, ApJ, 713, 837
  * 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, A147
  * Bond et al. (2001) Bond, I. A., Abe, F., Dodd, R. J., et al. 2001, MNRAS, 327, 868
  * Bessell & Brett (1988) Bessell, M. S., & Brett, J. M. 1988, PASP, 100, 1134
  * Danĕk & Heyrovský (2015a) Danĕk, K., & Heyrovský, D. 2015a, ApJ, 806, 63
  * Danĕk & Heyrovský (2015b) Danĕk, K., & Heyrovský, D. 2015b, 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
  * Dong et al. (2008) Dong, S., DePoy, D. L., Gaudi, B. S., et al. 2006, ApJ, 642, 842
  * Gaudi et al. (2008) Gaudi, B. S., Bennett, D. P., Udalski, A., et al. 2008, Science, 319, 927
  * Gaudi et al. (1998) Gaudi, B. S., Naber, R. M., & Sackett, P. D. 1998, ApJ, 502, L33
  * Gould (1992) Gould, A. 1992, ApJ, 392, 442
  * Gould (2000) Gould, A. 1992, ApJ, 542, 785
  * Gould et al. (2014) Gould, A., Udalski, A., Shin, I. -G., et al. 2014, Science, 345, 46
  * Griest & Hu (1993) Griest, K., & Hu, W. 1993, ApJ, 407, 440
  * Han et al. (2019) Han, C., Bennett, D. P., Udalski, A., et al. 2019, AJ, 158, 114
  * Han & Gould (1997) Han, C., & Gould, A. 1997, ApJ, 480, 196
  * Han et al. (2022a) Han, C., Gould, A., Bond, I. A., et al. 2022a, A&A, 662, A70
  * Han et al. (2022a) Han, C., Kim, D., Jung, Y. K. 2020a, AJ, 160, 17
  * Han et al. (2022b) Han, C., Kim, D., Yang, H. 2022b, A&A, 664, 114
  * Han et al. (2023) Han, C., Lee, C.-U., Gould, A., et al. 2023, A&A, in press
  * Han et al. (2020b) Han, C., Lee, C.-U., Udalski, A., et al. 2020b, AJ, 159, 48
  * Han et al. (2022c) Han, C., Ryu, Y.-H., Shin, I.-G., et al. 2022c, A&A, 667, A64
  * Han et al. (2013) Han, C., Udalski, A., Choi, J.-Y., et al. 2013, ApJ, 762, L28
  * Han et al. (2021a) Han, C., Udalski, A., Kim, D., et al. 2021a, AJ, 161, 270
  * Han et al. (2021b) Han, C., Udalski, A., Lee, C.-U. 2021b, AJ, 162, 203
  * Han et al. (2022d) Han, C., Udalski, A., Lee, C.-U. 2022d, A&A, 658, A93
  * Hwang et al. (2018) Hwang, K. -H., Udalski, A., Shvartzvald, Y., et al. 2018, AJ, 155, 20
  * Jung et al. (2021) Jung, Y. K., Han, C., Udalski, A., et al. 2021, AJ, 161, 293
  * Kervella et al. (2004) Kervella, P., Thévenin, F., Di Folco, E., & Ségransan, D. 2004, A&A, 426, 29
  * Kim et al. (2018a) Kim, D.-J., Kim, H.-W., Hwang, K.-H., et al. 2018a, AJ, 155, 76
  * Kim et al. (2018b) Kim, H.-W., Hwang, K.-H., Shvartzvald, Y., et al. 2018b, arXiv:1806.07545
  * Kim et al. (2016) Kim, S.-L., Lee, C.-U., Park, B.-G., et al. 2016, JKAS, 49, 37
  * Kuang et al. (2022) Kuang, R., Zang, W., Jung, Y. K., et al. 2022, MNRAS, 516, 1704
  * Madsen & Zhu (2019) Madsen, S., & Zhu, W. 2019, ApJ, 878, L29
  * Mao & Paczyński (1991) Mao, S., & Paczyński, B., 1991, ApJ, 374, L37
  * Nataf et al. (2013) Nataf, D. M., Gould, A., Fouqué, P., et al. 2013, ApJ, 769, 88
  * Paczyński (1986) Paczyński, B. 1986, ApJ, 301, 503
  * Poleski et al. (2014) Poleski, R., Skowron, J., Udalski, A., et al. 2014, ApJ, 795, 42
  * Ryu et al. (2020) Ryu, Y.-H., Udalski, A., Yee, J. C., et al. 2020, AJ, 160, 183
  * Schneider & Weiss (1986) Schneider, P., & Weiss, A. 1986, A&A, 164, 237
  * Schneider & Weiss (1987) Schneider, P., & Weiss, A. 1987, A&A, 171, 49
  * Suzuki et al. (2018) Suzuki, D., Bennett, D. P., Udalski, A., et al. 2018, AJ, 155, 263
  * Udalski et al. (1994) Udalski, A., Szymański, M., Kałużny, J., et al. 1994, Acta Astron., 44, 1
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1
  * 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
  * Zang et al. (2021) Zang, W., Han, C., Kondo, I., et al. 2021, Res. Astron. Astrophys., 21, 239