  Given a $C^{1,\beta}$ regular open set $D$ with boundary $\partial D =
\bigcup_{k=1}^d \bigcup_{j=1}^{m_k} \Gamma_{k,j}$, where for any $1\le k\le d$,
$ 1 \le j \le m_k$, $\Gamma_{k,j}$ is a $C^{1,\beta}$ submanifold without
boundary of codimension $k$ and $\{\Gamma_{k,j}\}_{1\le k\le d, 1\le j \le
m_k}$ are disjoint. We show that the heat kernel $p^D(t,x,y)$ of the following
non-local operator with multi-singular critical killing potential
\begin{align*} \big((\Delta|_D)^{\alpha/2} - \kappa\big)(f)(x):=
\text{p.v.}\mathcal{A}(d,-\alpha) \int_D \frac{f(y)-f(x)}{|y-x|^{d+\alpha}}dy -
\sum_{k=1}^d \sum_{j=1}^{m_k} \lambda_{k,j}
\delta_{\Gamma_{k,j}}(x)^{-\alpha}f(x), \end{align*} where $ \lambda_{k,j}>0,
\alpha \in (0,2), \beta \in ((\alpha-1)_+,1]$, has the following estimates: for
any given $T>0$, \begin{align*} p^D(t,x,y) \asymp p(t,x,y) \prod_{k=1}^d
\prod_{j=1}^{m_k} (\frac{\delta_{\Gamma_{k,j}}(x)}{t^{1/\alpha}} \wedge
1)^{p_{k,j}}(\frac{\delta_{\Gamma_{k,j}}(y)}{t^{1/\alpha}} \wedge 1)^{p_{k,j}},
\forall t\in (0,T), x,y\in D, \end{align*} where $p(t,x,y)$ is the heat kernel
of the $\alpha$-stable process on $\mathbb R^d$, and $p_{k,j}$ are determined
by $\lambda_{k,j}$ via a strictly increasing function $\lambda =
C(d,k,\alpha,p)$. Our method is based on the result established in [Cho et al.
Journal de Math\'ematiques Pures et Appliqu\'ees 143(2020): 208-256] and a
detailed analysis of $C^{1,\beta}$ manifolds.
