# Spectral convergence of the Laplace operator with Robin boundary conditions
on a small hole

Diana Barseghyan111corresponding author and Baruch Schneider

(Department of Mathematics, University of Ostrava, 30.dubna 22, Ostrava 70103,
Czech Republic
E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>)

###### Abstract

In this paper we study a bounded domain with a small hole removed. Our main
result concerns the spectrum of the Laplace operator with the Robin conditions
imposed at the hole boundary. Moreover we prove that under some suitable
assumptions on the parameter in the boundary condition the spectrum of the
Laplacian converges in the Hausdorff distance sense to the spectrum of the
Laplacian defined on the unperturbed domain.

Keywords. Robin Laplacian, spectral convergence, domain with a hole.

Mathematics Subject Classification. Primary: 58J50; Secondary: 35P15, 47A10.

##  1 Introduction

It is a common expectation that small perturbations of the physical situation
will lead to only a small change in the spectrum. In the case of domain
perturbations, this is largely true for Dirichlet boundary conditions, while
the Neumann or Robin case is more delicate. In the recent literature, such
questions have already received quite a few answers, starting with the seminal
work of Rauch and Taylor on the spectrum of the Laplace operator of domains
with holes [13].

An excellent shortcut to the recent work on the asymptotic behavior of the
eigenvalues of the Laplace operator on the domain with small spherical
obstacles imposing the Neumann condition at their boundary and the Dirichlet
condition at the rest part of the boundary can be found in [11].

Maz’ya, Nazarov and Plamenewskii, see [10, Ch.9, vol.I], have considered the
Laplace operator on the domain with obstacles, imposing the Dirichlet
condition on their boundary and have proved the validity of a complete
asymptotic expansion for the eigenvalues.

For a survey on more recent research in this subject, we refer the reader to
[4], [5], [6] where authors have considered the Dirichlet Laplacians on
Euclidean domains or manifolds with holes and studied the problems of the
resolvent convergence.

The problems with small Neumann obstacles of more general geometry can cause
abrupt changes in the spectrum. For example, such an effect is observed when
the hole has a ”split ring” geometry, see [14]. The split ring (even if very
small) can produce additional eigenvalues that have nothing in common with the
eigenvalues of the Neumann Laplacian on the unperturbed domain. The problems
with the Neumann obstacles having more general geometry have been studied in
[5] and later in [2].

The Robin case for general self-adjoint elliptic operators was considered in
[8] but with the restriction that the boundary of the unperturbed domain is
$C^{2}$\- smooth. In the mentioned work shrinking the hole and scaling
properly the parameter in the boundary condition, the authors obtain an
operator family that converges, in the norm-resolvent sense, to an operator
with a point interaction in the domain without the hole. Results on resolvent
convergence for operators with Robin conditions in domains with small holes in
higher dimensions were also considered in [3], but in this paper the original
domain must again be $C^{2}$-smooth and the Robin condition was used with the
coefficient independent of the hole size.

In this paper we will focus our attention on the spectral properties of a
Laplacian defined on a two-dimensional bounded domain with no additional
assumptions on the smoothness of its boundary with a single hole
$K_{\varepsilon}$ (for a fixed parameter $\varepsilon$) having the Lipschitz
boundary. On the boundary of the original domain we impose the Dirichlet or
Neumann condition, and on the boundary of the hole we impose the Robin
condition with the coefficient depending on the size of the hole.

Our main result is the proof of the spectral convergence, in the Hausdorff
distance sense, of the spectrum of the Laplacian defined on the perturbed
domain to the spectrum of the Laplacian defined on the original domain.

Plan of the paper. The paper consists of 6 sections, besides this
introduction.

In Section 2 we present the main results, and consider a general theorem,
namely Theorem 2.4. We will use Theorem 2.4 in the proof of Theorem 2.1 about
the spectral convergence for the Laplacians on $\Omega$ and $\Omega\setminus
K_{\varepsilon}$.

Section 3 contains the main tools of the spectral convergence of operators on
varying Hilbert spaces

In Section 4 and Section 5 we prove our results to which we already alluded.

In Section 6 we give some auxiliary material established in [2].

##  2 Main results

In this section, we present our main results. These results are proven in the
following sections.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain and $K\subset\Omega$ be
a compact simply connected set with Lipschitz boundary. We denote
$\Omega_{K}:=\Omega\setminus K$. By using Lemma 4.2, the quadratic forms

$\displaystyle q_{1}(u)=\int_{\Omega_{K}}|\nabla
u|^{2}\,dx\,dy+\gamma_{K}\int_{\partial K}|u|_{K}^{2}\,d\mu,\quad
u\in\mathcal{H}^{1}(\Omega_{K})\,,$ (2.1) $\displaystyle
q_{2}(u)=\int_{\Omega_{K}}|\nabla u|^{2}\,dx\,dy+\gamma_{K}\int_{\partial
K}|u|_{K}^{2}\,d\mu,\quad
u\in\mathcal{H}^{1}(\Omega_{K})\cap\mathcal{H}_{0}^{1}(\Omega)\,,$ (2.2)

where $\mu$ is the measure on $\partial K$ related to the arc length and
$\gamma_{K}$ is a real number, are closed and semi-bounded from below and
hence define unique self-adjoint operators $H_{\Omega_{K}}(\gamma_{K})$ and
$\widetilde{H}_{\Omega_{K}}(\gamma_{K})$ which act as the Laplacian on their
domains. We will study the question of the convergence of the spectrums of the
operators $H_{\Omega_{K}}(\gamma_{K}),\widetilde{H}_{\Omega_{K}}(\gamma_{K})$
when the hole $K$ converges to a point.

We start by a rather important results in the following theorem.

######  Theorem 2.1.

Let $\Omega$ be an open bounded domain in $\mathbb{R}^{2}$. Suppose that
${\mathbb{B}}_{\varepsilon}\subset\Omega$ is a ball with center at some point
$x_{0}\in\Omega$ and radius $\varepsilon>0$. Suppose that
$K=K_{\varepsilon}\subset\mathbb{B}_{\varepsilon}$ be a bounded simply
connected compact set with Lipschitz boundary. Let $H^{N}_{\Omega}$ and
$H^{D}_{\Omega}$ be the Neumann and Dirichlet Laplacians defined on the
unperturbed domain $\Omega$ and
$H_{\Omega_{K_{\varepsilon}}},\widetilde{H}_{\Omega_{K_{\varepsilon}}}$ be the
operators generated by (2.1) and (2.2) on $\Omega_{K_{\varepsilon}}$ with the
coefficient $\gamma_{K_{\varepsilon}}>0$ satisfying

$\gamma_{K_{\varepsilon}}=\varepsilon M_{\varepsilon},$ (2.3)

where $M_{\varepsilon}$ explodes to infinity under the condition that
$M_{\varepsilon}={o}\left(\frac{1}{\varepsilon^{3/2}}\right)$. Then, for
sufficiently small $\varepsilon$, there exists $\eta(\varepsilon)>0$ with
$\eta(\varepsilon)\to 0$ as $\varepsilon\to 0$, so that the following spectral
convergence occurs

$\overline{d}\left(\sigma_{\bullet}\left(H_{\Omega_{K_{\varepsilon}}}\right),\,\sigma_{\bullet}\left(H^{N}_{\Omega}\right)\right)\leq\eta(\varepsilon),$
$\overline{d}\left(\sigma_{\bullet}\left(\widetilde{H}_{\Omega_{K_{\varepsilon}}}\right),\,\sigma_{\bullet}\left(H^{D}_{\Omega}\right)\right)\leq\eta(\varepsilon),$

where $\overline{d}$ is defined in (3.9) and $\sigma_{\bullet}(\cdot)$ denotes
either the entire spectrum, the essential spectrum, or the discrete spectrum.
Furthermore, the multiplicity of the discrete spectrum is preserved.

The previous result motivates the following consequences:

######  Corollary 2.2.

Suppose that $H^{N}_{\Omega}$ has purely discrete spectrum denoted by
$\lambda^{N}_{k}(\Omega)$ (repeated according to multiplicity), and let
$\lambda^{D}_{k}(\Omega)$ be the discrete spectrum of $H^{D}_{\Omega}$. Then
the infimum of the essential spectrums of
$H_{\Omega_{K_{\varepsilon}}},\widetilde{H}_{\Omega_{K_{\varepsilon}}}$ tend
to infinity and there exists $\eta_{k}(\varepsilon)>0$ with
$\eta_{k}(\varepsilon)\to 0$ as $\varepsilon\to 0$ such that

$|\lambda^{N}_{k}(\Omega)-\lambda_{k}(\Omega_{K_{\varepsilon}})|\leq\eta_{k}(\varepsilon),$
$|\lambda^{D}_{k}(\Omega)-\beta_{k}(\Omega_{K_{\varepsilon}})|\leq\eta_{k}(\varepsilon)$

for small enough $\varepsilon$. Here, $\lambda_{k}(\Omega_{K_{\varepsilon}})$
and $\beta_{k}(\Omega_{K_{\varepsilon}})$ denote the discrete spectrum of
$H_{\Omega_{K_{\varepsilon}}}$ and $\widetilde{H}_{\Omega_{K_{\varepsilon}}}$
(below the essential spectrum) repeated according to multiplicity.

######  Corollary 2.3.

The Hausdorff distance between the spectra of $H_{\Omega_{K_{\varepsilon}}}$
and $H^{N}_{\Omega}$ and $\widetilde{H}_{\Omega_{K_{\varepsilon}}}$ and
$H^{D}_{\Omega}$ converges to zero on any compact interval $[0,\Lambda]$.

The proof of Theorem 2.1 is based on Theorem 3.1 and on the following theorem:

######  Theorem 2.4.

Under the assumptions of Theorem 2.1 the operators $H^{N}_{\Omega}$ and
$H_{\Omega_{K_{\varepsilon}}}$ are $\delta(\varepsilon)$ close of order $2$
with $\delta(\varepsilon)\to 0$ as $\varepsilon\to 0$. The same is true for
the operators $H^{D}_{\Omega}$ and $\widetilde{H}_{\Omega_{K_{\varepsilon}}}$.

##  3 Main tool of the spectral convergence of operators on varying Hilbert
spaces

For the convenience of the reader, this section begins by reviewing some basic
facts that ensure spectral convergence for two operators having different
domains. For more information we refer the reader to [12]. To a Hilbert space
$H$ with inner product $(\cdot,\cdot)$ and norm $\|\cdot\|$ together with a
non-negative, unbounded operator $A$ we associate the scale of Hilbert spaces

$H_{k}:=\mathrm{Dom}((A+I)^{k/2}),\quad\|u\|_{k}:=\|(A+I)^{k/2}u\|,\,\,k\geq
0,$

where $I$ is the identity operator.

We think of $(H^{\prime},A^{\prime})$ as some perturbation of $(H,A)$ and want
to relax the assumption so that the spectral properties are not the same, but
still close.

######  Definition 1.

(see [12]) Suppose we have linear operators

$\displaystyle J:H\longrightarrow H^{\prime},\quad\quad
J_{1}:H_{1}\longrightarrow H_{1}^{\prime}$ $\displaystyle
J^{\prime}:H^{\prime}\longrightarrow H,\quad\quad
J_{1}^{\prime}:H_{1}^{\prime}\longrightarrow H_{1}.$

Let $\delta>0$ and $k\geq 1$. We say that $(H,A)$ and
$(H^{\prime},A^{\prime})$ are $\delta$-close of order $k$ iff the following
conditions are fulfilled:

$\displaystyle\|Jf-J_{1}f\|_{0}\leq\delta\|f\|_{1},$ (3.1)
$\displaystyle|(Jf,u)-(f,J^{\prime}u)|\leq\delta\|f\|_{0}\|u\|_{0},$ (3.2)
$\displaystyle\|u-JJ^{\prime}u\|_{0}\leq\delta\|u\|_{1},$ (3.3)
$\displaystyle\quad\|Jf\|_{0}\leq 2\|f\|_{0},\,\,\,\|J^{\prime}u\|_{0}\leq
2\|u\|_{0},$ (3.4) $\displaystyle\|(f-J^{\prime}Jf)\|_{0}\leq\delta\|f\|_{1},$
(3.5) $\displaystyle\|J^{\prime}u-J_{1}^{\prime}u\|_{0}\leq\delta\|u\|_{1}$
(3.6)
$\displaystyle|a(f,J_{1}^{\prime}u)-a^{\prime}(J_{1}f,u)|\leq\delta\|f\|_{k}\|u\|_{1},$
(3.7)

for all $f,u$ in the appropriate spaces. Here, $a$ and $a^{\prime}$ denote the
sesquilinear forms associated to $A$ and $A^{\prime}$.

We denote by $d_{\mathrm{Haussdorff}}(A,B)$ the Hausdorff distance for subsets
$A,B\subset\mathbb{R}$:

$d_{\mathrm{Haussdorff}}(A,B):=\mathrm{max}\left\\{\underset{a\in
A}{\mathrm{sup}}\,d(a,B),\,\underset{b\in B}{\mathrm{sup}}\,d(b,A)\right\\},$
(3.8)

where $d(a,B):=\mathrm{inf}_{b\in B}|a-b|$. We set

$\overline{d}(A,B):=d_{\mathrm{Hausdorff}}\left((A+I)^{-1},(B+I)^{-1}\right)$
(3.9)

for closed subsets of $[0,\infty)$.

For the next result, which comes from the work of O. Post [12], we have the
following spectral convergence theorem in terms of the distance
$\overline{d}$.

######  Theorem 3.1.

[12] There exists $\eta(\delta)>0$ with $\eta(\delta)\to 0$ as $\delta\to 0$
such that

$\overline{d}(\sigma_{\bullet}(A),\,\sigma_{\bullet}(A^{\prime}))\leq\eta(\delta)$
(3.10)

for all pairs of non-negative operators and Hilbert spaces $(H,A)$ and
$(H^{\prime},A^{\prime})$ that are $\delta$-close. Here, $\sigma_{\bullet}(A)$
denotes either the entire spectrum, the essential or the discrete spectrum of
$A$. Also, the multiplicity of the discrete spectrum,
$\sigma_{\mathrm{disc}}$, is preserved, i.e. if
$\lambda\in\sigma_{\mathrm{disc}}(A)$ has multiplicity $\mu>0$, then there
exist $\mu$ eigenvalues (not necessarily all distinct) of the operator
$A^{\prime}$ belonging to the interval
$(\lambda-\eta(\delta),\lambda+\eta(\delta))$.

We now turn to the proof of Theorem 2.4. Since the proof is almost the same
for both the Dirichlet and Neumann cases, we will restrict ourselves to the
Neumann case. The only difference is Lemma 6.3, but its validity for the
Dirichlet case can be easily checked from [2].

##  4 Proof of Theorem 2.4

At this stage, we divide the proof into two steps.

Step 1. Construction of the mappings $J,J^{\prime},J_{1},J_{1}^{\prime}$.

We can apply the technique of [2]. It is easy to see that
$H=L^{2}(\Omega),\,H^{\prime}=L^{2}(\Omega_{K_{\varepsilon}}),\,A=A^{\prime}=-\Delta$,
$H_{1},\,H_{1}^{\prime}$ correspond to Sobolev spaces
$\mathcal{H}^{1}(\Omega)$ and $\mathcal{H}^{1}(\Omega_{K_{\varepsilon}})$ and
$H_{2}=\mathrm{Dom}(H^{N}_{\Omega})$. The norm $\|\cdot\|_{0}$ corresponds
with the $L^{2}$ norm,

$\displaystyle\|u\|_{1}=\left(\|u\|^{2}_{0}+\|\nabla
u\|_{0}^{2}+\gamma_{K_{\varepsilon}}\int_{\partial
K_{\varepsilon}}|u|^{2}\,d\mu\right)^{1/2},\quad\text{for}\quad
u\in\mathcal{H}^{1}(\Omega_{K_{\varepsilon}}),$
$\displaystyle\|u\|_{1}=\left(\|u\|^{2}_{0}+\|\nabla
u\|_{0}^{2}\right)^{1/2},\quad\text{for}\quad u\in\mathcal{H}^{1}(\Omega)$
$\displaystyle\text{and}\quad\|f\|_{2}=\|-\Delta f+f\|_{0}.$

We define $Ju=J_{1}u=u|_{\chi_{\Omega_{K_{\varepsilon}}}}$ for all $u\in H$
and $J^{\prime}u=u\chi_{\Omega_{K_{\varepsilon}}}$ for all $u\in H^{\prime}$.

Now let us construct the mapping $J_{1}^{\prime}:H_{1}^{\prime}\to H_{1}$.
Without loss of generality, assume that the ball $\mathbb{B}_{\varepsilon}$
mentioned in Theorem 2.1 and Theorem 2.4 is centered at the origin. Let
$\epsilon\in(\varepsilon,2\varepsilon)$ be a number to be chosen later and let
$\mathbb{B}_{\epsilon}\supset\mathbb{B}_{\varepsilon}$ be the ball again
centered on the origin and radius $\epsilon$,
$\Omega_{\epsilon}:=\Omega\backslash\mathbb{B}_{\epsilon}$.

We will first construct the mapping $J_{1}^{\prime}$ first for smooth
functions. For each $v\in C^{\infty}(\Omega_{K_{\varepsilon}})$ we define

$J_{1}^{\prime}v:=\begin{cases}v,\quad\text{on}\quad\Omega_{\epsilon},\\\
\frac{r}{\epsilon}\tilde{v}(\epsilon,\varphi),\quad\text{on}\quad\mathbb{B}_{\epsilon},\end{cases}$

where $\tilde{v}(r,\varphi)=v(r\cos\varphi,r\sin\varphi)$.

Now let us construct the mapping $J_{1}^{\prime}u$ for any $u\in
H_{1}^{\prime}$. Using the approximation method described in [7, Thm.2,
5.3.2], for the fixed sequence $\\{\eta_{k}\\}_{k=1}^{\infty}$ converging to
zero we construct the sequence $v_{\eta_{k}}\in
C^{\infty}(\Omega_{K_{\varepsilon}})$ which satisfies

$\int_{\Omega_{K_{\varepsilon}}}|\nabla(u-v_{\eta_{k}})|^{2}\,dxdy+\int_{\Omega_{K_{\varepsilon}}}|u-v_{\eta_{k}}|^{2}\,dxdy<\eta_{k}\,\|u\|_{1}^{2}.$
(4.11)

To deal with $\int_{\partial K_{\varepsilon}}|u|^{2}\,d\mu$ we will use the
trace inequality [7]. We present it immediately after an auxiliary result on
Lipschitz bounds [7].

######  Lemma 4.1.

Let $\omega$ be a bounded two-dimensional open set with a Lipschitz boundary.
Then there exist $\delta>0$ and $\mu\in C^{\infty}(\overline{\omega})$ such
that

$\mu\cdot\nu\geq\delta,\quad\text{a.e. on}\quad\partial\omega,$

where $\nu$ is the normal to $\partial\omega$.

######  Lemma 4.2.

Let $v\in\mathcal{H}^{1}(\omega)$, where $\omega$ is a two-dimensional domain
with Lipschitz boundary. Then there exists a constant $K>0$ depending on
$\omega$ such that

$\int_{\partial\omega}|v|^{2}d\mu\leq K\int_{\omega}\left(\delta|\nabla
v|^{2}+\frac{1}{\delta}|v|^{2}\right)\,dx\,dy,$

where $K$ depends only on the norm of $\mu$ in $C^{1}(\overline{\Omega})$ and
$\delta\in(0,1)$.

Combining the above lemma with $\delta=\frac{1}{2}$ and the inequality (4.11)
we get

$\displaystyle\int_{\partial K_{\varepsilon}}|u-v_{\eta_{k}}|^{2}\,d\mu$
$\displaystyle\leq
2K\,\left(\int_{\Omega_{K_{\varepsilon}}}|\nabla(u-v_{\eta_{k}})|^{2}\,dx\,dy+\int_{\Omega_{K_{\varepsilon}}}|u-v_{\eta_{k}}|^{2}\,dx\,dy\right)$
$\displaystyle<2K\eta_{k}\|u\|_{1}^{2}.$

Therefore

$\|u-v_{\eta_{k}}\|_{1}<\sqrt{(1+2K\gamma_{K_{\varepsilon}})\,\eta_{k}}\,\|u\|_{1}.$
(4.12)

Let us mention that in view of the inequalities (4.19) and (4.25), which will
be proved later, and the construction of the function $J_{1}^{\prime}$, it
follows that for any smooth function $v$ the integrals $\int_{\Omega}|\nabla
J_{1}^{\prime}v|^{2}\,dx\,dy$ and $\int_{\Omega}|J_{1}^{\prime}v|^{2}dx\,dy$
can be estimated from above by $\|v\|_{1}^{2}$ multiplied by some constant.
Combining this with the Lemma 4.2 we get

$\|J_{1}^{\prime}v\|_{1}^{2}\leq\overline{C}(\varepsilon)\|v\|_{1}^{2},$
(4.13)

where $\overline{C}(\varepsilon)$ is some constant.

Due to the positivity of the coefficient $\gamma_{K_{\varepsilon}}$ and Lemma
4.2 the completeness of the space $H_{1}$ is equivalent to the completeness of
the Sobolev space $\mathcal{H}^{1}$. Thus, using (4.12) and (4.13) we can
define

$J_{1}^{\prime}u=\underset{k\to\infty}{\mathrm{lim}}\,J_{1}^{\prime}v_{\eta_{k}}.$
(4.14)

Step 2. The conditions (3.1)-(3.7) hold for the mappings
$J,J^{\prime},J_{1},J_{1}^{\prime}$.

Indeed, we have that the estimates (3.1)-(3.4) are satisfied with $\delta=0$.

We now prove (3.5), i.e., under the assumptions stated in Theorem 2.4,
inequality (3.5) is satisfied with
$\delta=\mathcal{O}\left(\sqrt{\frac{\varepsilon}{\gamma_{K_{\varepsilon}}}+\varepsilon}\right)$
for small enough $\varepsilon$. Thus, given (2.3), $\delta$ converges to zero
as $\varepsilon\to 0$.

Given our construction, we have

$\displaystyle\|f-J^{\prime}Jf\|_{0}^{2}=\int_{\Omega}|f-J^{\prime}f|_{\Omega_{K_{\varepsilon}}}|^{2}\,dx\,dy=\int_{\Omega}|f-f\chi_{\Omega_{K_{\varepsilon}}}|^{2}\,dx\,dy$
$\displaystyle=\int_{K_{\varepsilon}}|f|^{2}\,dxdy\leq\int_{\mathbb{B}_{\varepsilon}}|f|^{2}\,dx\,dy.$
(4.15)

To complete the proof of (3.5), we use the following lemma, applied with
$\eta=\varepsilon$ and $\Gamma=\emptyset$:

######  Lemma 4.3.

Let $\Omega$ be an open bounded domain in $\mathbb{R}^{2}$. Suppose that
$\mathbb{B}_{\epsilon}\subset\Omega$ is a ball with center at some point
$x_{0}\in\Omega$ and radius $\epsilon>0$. Suppose that
$\Gamma\subset\mathbb{B}_{\epsilon}$ be a bounded simply connected compact set
with Lipschitz boundary. Then for any function
$u\in\mathcal{H}^{1}(\mathbb{B}_{\epsilon}\setminus\Gamma)$ the following
inequality holds

$\int_{\mathbb{B}_{\epsilon}\setminus\Gamma}|u(x,y)|^{2}\,dxdy\leq
C_{1}\left(\frac{\epsilon}{\gamma_{\Gamma}}+\epsilon\right)\|u\|_{1}^{2},$

where $C_{1}>0$ is a constant that depends on the distance between the
boundary of $\mathbb{B}_{\epsilon}$ and the boundary of $\Omega$.

Let us pass to (3.6), i.e. under the assumptions given in Theorem 2.4
inequality (3.6) is satisfied with
$\delta=\mathcal{O}\left(\sqrt{\frac{\varepsilon}{\gamma_{K_{\varepsilon}}}+\varepsilon}\right)$
for small enough $\varepsilon$. Thus, in view of (2.3), $\delta$ converges to
zero as $\varepsilon\to 0$.

Using the construction of $J^{\prime}$ and $J_{1}^{\prime}$ and the
completeness of $C^{\infty}(\Omega_{K_{\varepsilon}})$ in space
$H_{1}^{\prime}$ it is sufficient to prove (3.6) for $u\in
C^{\infty}(\Omega_{K_{\varepsilon}})$.

Considering that $J^{\prime}u=0$ on $K_{\varepsilon}$, $J^{\prime}u=u$ on
$\Omega_{K_{\varepsilon}}$ and $J_{1}^{\prime}u=u$ on
$\Omega\setminus\mathbb{B}_{\epsilon}$ one has

$\displaystyle\|J^{\prime}u-J_{1}^{\prime}u\|_{0}^{2}=\int_{\Omega_{K_{\varepsilon}}}|u-J_{1}^{\prime}u|^{2}\,dxdy+\int_{K_{\varepsilon}}|J_{1}^{\prime}u|^{2}\,dxdy$
$\displaystyle=\int_{\mathbb{B}_{\epsilon}\setminus
K_{\varepsilon}}|u-J_{1}^{\prime}u|^{2}\,dxdy+\int_{K_{\varepsilon}}|J_{1}^{\prime}u|^{2}\,dxdy$
$\displaystyle\leq 2\int_{\mathbb{B}_{\epsilon}\setminus
K_{\varepsilon}}|u|^{2}\,dxdy+2\int_{\mathbb{B}_{\epsilon}}|J_{1}^{\prime}u|^{2}\,dxdy+\int_{K_{\epsilon}}|J_{1}^{\prime}u|^{2}\,dxdy.$
(4.16)

To estimate the first integral on the right-hand side of (4.16), we use the
Lemma 4.3 with $\eta=\epsilon$ and $\Gamma=K_{\varepsilon}$.

Now let us examine the second term. Passing to polar coordinates, we get

$\displaystyle\int_{\mathbb{B}_{\epsilon}}|J_{1}^{\prime}u|^{2}\,dx\,dy=\int_{0}^{\epsilon}\int_{0}^{2\pi}r\left|\frac{r}{\epsilon}\tilde{u}(\epsilon,\varphi)\right|^{2}\,dr\,d\varphi$
$\displaystyle\leq\int_{0}^{\epsilon}\int_{0}^{2\pi}r|\tilde{u}(\epsilon,\varphi)|^{2}\,dr\,d\varphi\leq\epsilon\int_{0}^{\epsilon}\int_{0}^{2\pi}|\tilde{u}(\epsilon,\varphi)|^{2}\,dr\,d\varphi$
$\displaystyle=\epsilon^{2}\int_{0}^{2\pi}|\tilde{u}(\epsilon,\varphi)|^{2}\,d\varphi.$

Taking into account that
$\epsilon\int_{0}^{2\pi}|\tilde{u}(\epsilon,\varphi)|^{2}\,d\varphi$ coincides
with the curvilinear integral
$\int_{\partial\mathbb{B}_{\epsilon}}|u|^{2}\,d\mu$, the above bound leads to

$\int_{\mathbb{B}_{\epsilon}}|J_{1}^{\prime}u|^{2}\,dxdy\leq\epsilon\int_{\partial\mathbb{B}_{\epsilon}}|u|^{2}\,d\mu.$
(4.17)

In view of Lemma 4.2 applied with $\delta=\frac{1}{2}$, one estimates from
above the right-hand side of (4.17) by
$2K\epsilon\int_{\Omega\setminus\mathbb{B}_{\epsilon}}(|\nabla
u|^{2}+|u|^{2})\,dx\,dy$. Then use the following obvious bound which holds due
to the positivity of $\gamma_{K_{\varepsilon}}$:

$\int_{\Omega\setminus\mathbb{B}_{\epsilon}}(|\nabla
g|^{2}+|g|^{2})\,dx\,dy\leq\|g\|_{1}^{2},\quad\text{for all}\quad
g\in\mathcal{H}^{1}(\Omega_{K_{\epsilon}}),$ (4.18)

and the fact that $\epsilon\leq 2\varepsilon$ we have

$\int_{\mathbb{B}_{\epsilon}}|J_{1}^{\prime}u|^{2}\,dxdy\leq
4K\varepsilon\|u\|_{1}^{2}.$ (4.19)

Indeed, we have

$\int_{K_{\varepsilon}}|J_{1}^{\prime}u|^{2}\,dx\,dy\leq
4K\varepsilon\|u\|_{1}^{2},$

since

$\int_{K_{\varepsilon}}|J_{1}^{\prime}u|^{2}\,dx\,dy\leq\int_{B_{\epsilon}}|J_{1}^{\prime}u|^{2}\,dx\,dy.$

Combining the above inequality together with (4.19), the right-hand side of
(4.16) can be estimated as follows

$\|J^{\prime}u-J_{1}^{\prime}u\|_{0}^{2}\leq\left(2C_{1}\left(\frac{\varepsilon}{\gamma_{K_{\varepsilon}}}+\varepsilon\right)+12K\varepsilon\right)\|u\|_{1}^{2}$

which proves (3.6) with
$\delta=\mathcal{O}\left(\sqrt{\frac{\varepsilon}{\gamma_{K_{\varepsilon}}}+\varepsilon}\right)$
for small enough $\varepsilon$.

We now give the proof of the estimate (3.7), i.e. under the assumptions given
in Theorem 2.4, the inequality (3.7) holds with $k=2$ and
$\delta=\mathcal{O}\left(\varepsilon^{1/6}+\gamma_{K_{\varepsilon}}^{1/2}\varepsilon^{1/4}\right)$
for sufficiently small $\varepsilon$. Thus, in view of (2.3), $\delta$
converges to zero as $\varepsilon\to 0$.

As before without loss of generality suppose that $u\in
C^{\infty}(\Omega_{K_{\varepsilon}})$. Since $J_{1}^{\prime}u=u$ on
$\Omega_{\mathbb{B}_{\epsilon}}$ and $J_{1}u=u$ on $\Omega_{K_{\varepsilon}}$
we have

$\displaystyle|a(f,J_{1}^{\prime}u)-a^{\prime}(J_{1}f,u)|$
$\displaystyle=\left|\int_{\Omega}\nabla
f\,\overline{\nabla(J_{1}^{\prime}u)}\,dx\,dy-\int_{\Omega_{K_{\varepsilon}}}\nabla(J_{1}f)\,\overline{\nabla
u}\,dx\,dy-\gamma_{K_{\varepsilon}}\int_{\partial
K_{\varepsilon}}(J_{1}f)\,\overline{u}\,d\mu\right|$
$\displaystyle=\biggl{|}\int_{\Omega_{\mathbb{B}_{\epsilon}}}\nabla
f\,\overline{\nabla u}\,dx\,dy+\int_{\mathbb{B}_{\epsilon}}\nabla
f\,\overline{\nabla(J_{1}^{\prime}u)}\,dx\,dy-\int_{\Omega_{K_{\varepsilon}}}\nabla(J_{1}f)\,\overline{\nabla
u}\,dx\,dy-\gamma_{K_{\varepsilon}}\int_{\partial
K_{\varepsilon}}(J_{1}f)\,\overline{u}\,d\mu\biggr{|}$
$\displaystyle=\biggl{|}\int_{\mathbb{B}_{\epsilon}}\nabla
f\,\overline{\nabla(J_{1}^{\prime}u)}\,dx\,dy-\int_{\mathbb{B}_{\epsilon}\backslash
K_{\varepsilon}}\nabla f\,\overline{\nabla
u}\,dx\,dy-\gamma_{K_{\varepsilon}}\int_{\partial
K_{\varepsilon}}(J_{1}f)\,\overline{u}\,d\mu\biggr{|}$
$\displaystyle\leq\left|\int_{\mathbb{B}_{\epsilon}}\nabla
f\,\overline{\nabla(J_{1}^{\prime}u)}\,dx\,dy\right|+\left|\int_{\mathbb{B}_{\epsilon}\backslash
K_{\varepsilon}}\nabla f\,\overline{\nabla
u}\,dx\,dy\right|+\gamma_{K_{\varepsilon}}\left|\int_{\partial
K_{\varepsilon}}(J_{1}f)\,\overline{u}\,d\mu\right|.$ (4.20)

Let us estimate each term of (4.20). Starting with the first term we get

$\left|\int_{\mathbb{B}_{\epsilon}}\nabla
f\,\overline{\nabla(J_{1}^{\prime}u)}\,dxdy\right|\leq\left(\int_{\mathbb{B}_{\epsilon}}|\nabla
f|^{2}\,dxdy\right)^{1/2}\left(\int_{\mathbb{B}_{\epsilon}}|\nabla(J_{1}^{\prime}u)|^{2}\,dxdy\right)^{1/2}.$
(4.21)

Since $f\in\mathcal{H}^{2}_{\mathrm{loc}}(\Omega)$, then using Lemma 6.1 (see
Appendix) applied with domain $\Omega^{\prime}$ such that
$\overline{\Omega^{\prime}}\subset\Omega$ and
$\mathbb{B}_{\epsilon}\subset\Omega^{\prime}$, and the fact that $\epsilon\leq
2\varepsilon$, the first term on the right hand side of (4.21) can be
estimated as follows

$\int_{\mathbb{B}_{\epsilon}}|\nabla f|^{2}\,dx\,dy\leq
2^{4/3}C_{2}\varepsilon^{4/3}\int_{\Omega^{\prime}}|-\Delta
f+f|^{2}\,dx\,dy\leq 2^{4/3}C_{2}\varepsilon^{4/3}\|f\|_{2}^{2}.$ (4.22)

To proceed further with the proof of an upper bound of (4.21), we need to
estimate the integral
$\int_{\mathbb{B}_{\epsilon}}|\nabla(J_{1}^{\prime}u)|^{2}\,dx\,dy$. Passing
to polar coordinates we get

$\displaystyle\int_{\mathbb{B}_{\epsilon}}|\nabla
J_{1}^{\prime}u|^{2}\,dx\,dy=\int_{\mathbb{B}_{\epsilon}}\left(\left|\frac{\partial(J_{1}^{\prime}u)}{\partial
x}\right|^{2}+\left|\frac{\partial(J_{1}^{\prime}u)}{\partial
y}\right|^{2}\right)\,dx\,dy$
$\displaystyle=\int_{0}^{\epsilon}\int_{0}^{2\pi}r\left|\frac{1}{\epsilon}\tilde{u}(\epsilon,\varphi)\,\cos\varphi-\frac{1}{\epsilon}\frac{\partial\tilde{u}}{\partial\varphi}(\epsilon,\varphi)\,\sin\varphi\right|^{2}\,dr\,d\varphi$
$\displaystyle+\int_{0}^{\varepsilon}\int_{0}^{2\pi}r\left|\frac{1}{\epsilon}\tilde{u}(\epsilon,\varphi)\,\sin\varphi+\frac{1}{\epsilon}\frac{\partial\tilde{u}}{\partial\varphi}(\epsilon,\varphi)\,\cos\varphi\right|^{2}\,dr\,d\varphi$
$\displaystyle\leq\frac{4}{\epsilon}\int_{0}^{\epsilon}\int_{0}^{2\pi}\left(|\tilde{u}(\epsilon,\varphi)|^{2}+\left|\frac{\partial\tilde{u}}{\partial\varphi}(\epsilon,\varphi)\right|^{2}\right)\,dr\,d\varphi$
$\displaystyle=4\int_{0}^{2\pi}|\tilde{u}(\epsilon,\varphi)|^{2}d\varphi+4\int_{0}^{2\pi}\left|\frac{\partial\tilde{u}}{\partial\varphi}(\epsilon,\varphi)\right|^{2}\,d\varphi.$
(4.23)

As in the proof of (3.6), we find that

$\int_{0}^{2\pi}|\tilde{u}(\epsilon,\varphi)|^{2}\,d\varphi=\frac{1}{\epsilon}\,\int_{\partial\mathbb{B}_{\epsilon}}|u|^{2}\,d\mu.$

Thus, using Lemma 4.2 applied with $\delta=\frac{1}{2}$ and inequality (4.18),
we have

$\int_{0}^{2\pi}|\tilde{u}(\epsilon,\varphi)|^{2}\,d\varphi\leq\frac{2K}{\epsilon}\|u\|_{1}^{2}.$
(4.24)

Now we come to the second term of (4.23). Given the Lemma 6.2 (see Appendix)
used for $g=u$, there exists a number $\tau\in(\varepsilon,2\varepsilon)$ such
that

$\int_{0}^{2\pi}\left|\frac{\partial\tilde{u}}{\partial\varphi}(\tau,\varphi)\right|^{2}\,d\varphi\leq
4\int_{\mathbb{B}_{2\varepsilon}\setminus\mathbb{B}_{\varepsilon}}(|\nabla
u|^{2}+|u|^{2})\,dxdy,$

where $\tilde{u}(r,\varphi):=u(r\cos\varphi,r\sin\varphi)$ and
$\tau\in(\varepsilon,2\varepsilon)$ is some number.

If $\tau$ belongs to the interval $(\varepsilon,3\varepsilon/2]$, then we take
$\epsilon$ as the supremum of all such numbers in
$(\varepsilon,3\varepsilon/2]$. In the opposite case if
$\tau\in(3\varepsilon/2,2\varepsilon)$, then let $\epsilon$ be the infimum of
such numbers. Since $u$ is a smooth function, the above inequality is
satisfied with $\tau=\epsilon$.

Combining this together with (4.18) we get

$\int_{0}^{2\pi}\left|\frac{\partial\tilde{u}}{\partial\varphi}(\epsilon,\varphi)\right|^{2}\,d\varphi\leq
4\|u\|_{1}^{2}.$

Thus, by virtue of (4.23), (4.24), the above estimate and using the fact that
$\epsilon\geq\varepsilon$, we have

$\int_{\mathbb{B}_{\epsilon}}|\nabla
J_{1}^{\prime}u|^{2}\,dx\,dy\leq\left(\frac{8K}{\varepsilon}+16\right)\|u\|_{1}^{2}.$
(4.25)

Finally, using the above bound and inequality (4.22), we estimate the right-
hand side of (4.21) as follows

$\left|\int_{\mathbb{B}_{\epsilon}}\nabla
f\,\overline{\nabla(J_{1}^{\prime}u)}\,dx\,dy\right|\leq
2^{13/6}(C_{2}(K+2\varepsilon))^{1/2}\varepsilon^{1/6}\|f\|_{2}\|u\|_{1}.$
(4.26)

Let us now consider the second term in (4.20). By virtue of (4.18) and (4.22)
we get

$\displaystyle\left|\int_{\mathbb{B}_{\epsilon}\backslash
K_{\varepsilon}}\nabla f\,\overline{\nabla
u}\,dx\,dy\right|\leq\left(\int_{\mathbb{B}_{\epsilon}}|\nabla
f|^{2}\,dx\,dy\right)^{1/2}\left(\int_{\mathbb{B}_{\epsilon}\backslash
K_{\varepsilon}}|\nabla u|^{2}\,dx\,dy\right)^{1/2}$ $\displaystyle\leq
2^{2/3}(C_{2})^{1/2}\varepsilon^{2/3}\|f\|_{2}\|u\|_{1}.$ (4.27)

Finally we move on to the third term in (4.20). We have

$\displaystyle\left|\int_{\partial
K_{\varepsilon}}(J_{1}f)\,\overline{u}\,d\mu\right|\leq\left(\int_{\partial
K_{\varepsilon}}|J_{1}f|^{2}\,d\mu\right)^{1/2}\left(\int_{\partial
K_{\varepsilon}}|u|^{2}\,d\mu\right)^{1/2}$
$\displaystyle=\left(\int_{\partial
K_{\varepsilon}}|f|^{2}\,d\mu\right)^{1/2}\left(\int_{\partial
K_{\varepsilon}}|u|^{2}\,d\mu\right)^{1/2}.$ (4.28)

Let us first find the appropriate estimate for the first integral on in the
right hand side of (4.28).

Let $\Pi(d)=(-d,d)^{2},\,d>0$, be the maximum square belonging to $\Omega$ and
containing $K_{\varepsilon}$. For almost all $x_{0}$ belonging to the
projection of $K_{\varepsilon}$ on the axis $X$, let $y(x_{0})\in(-d,d)$ be
the point such that

$|f(x_{0},y(x_{0}))|^{2}\leq\frac{1}{d}\int_{-d}^{d}|f(x_{0},t)|^{2}\,dt\,.$

Let us fix any $y\in(-d,d)$. Without loss of generality suppose that
$y>y(x_{0})$. Then

$\displaystyle|f(x_{0},y)|^{2}=\left|f(x_{0},y(x_{0}))+\int_{y(x_{0})}^{y}\frac{\partial
f}{\partial t}(x_{0},t)\,dt\right|^{2}$ $\displaystyle\leq
2\left|f(x_{0},y(x_{0}))\right|^{2}+2\left|\int_{y(x_{0})}^{y}\frac{\partial
f}{\partial t}(x_{0},t)\,dt\right|^{2}$
$\displaystyle\leq\frac{2}{d}\int_{-d}^{d}|f(x_{0},t)|^{2}\,dt+2(y-y(x_{0}))\int_{y_{0}(x_{0})}^{y}\left|\frac{\partial
f}{\partial t}(x_{0},t)\right|^{2}\,dt$
$\displaystyle\leq\frac{2}{d}\int_{-d}^{d}|f(x_{0},t)|^{2}\,dt+4d\int_{-d}^{d}\left|\frac{\partial
f}{\partial t}(x_{0},t)\right|^{2}\,dt\,.$ (4.29)

Without loss of generality assume that the boundary of the unperturbed set
which is $\frac{1}{\varepsilon}K_{\varepsilon}$ is parametrized as
$(x,y_{1}(x)),\,x\in(-1,1)$, where $y_{1}$ is some $C^{1}$-smooth function.
Then the parameterization of the boundary of $K_{\varepsilon}$ coincides with
$(x,\varepsilon y_{1}(x/\varepsilon)),\,x\in(-\varepsilon,\varepsilon)$.

Integrating $|f(x,y)|^{2}$ over $\partial K_{\varepsilon}$ and using the
inequality (4.29) we get

$\displaystyle\int_{\partial
K_{\varepsilon}}|f(x,y)|^{2}\,d\mu=\int_{-\varepsilon}^{\varepsilon}|f(x,\varepsilon
y_{1}(x/\varepsilon))|^{2}\,\left(1+y_{1}^{\prime
2}(x/\varepsilon)\right)^{1/2}\,dx$
$\displaystyle\leq\left(1+\|y_{1}^{\prime}\|_{L^{\infty}(-1,1)}^{2}\right)^{1/2}\int_{-\varepsilon}^{\varepsilon}\biggl{(}\frac{2}{d}\int_{-d}^{d}|f(x,t)|^{2}\,dt+4d\int_{-d}^{d}\left|\frac{\partial
f}{\partial t}(x,t)\right|^{2}\,dt\biggr{)}\,dx$ $\displaystyle\leq
2\mathrm{max}\left\\{1/d,2d\right\\}\left(1+\|y_{1}^{\prime}\|_{L^{\infty}(-1,1)}^{2}\right)^{1/2}\int_{\Omega_{\varepsilon}}(|f|^{2}+|\nabla
f|^{2})\,dxdy\,,$ (4.30)

where $\Omega_{\varepsilon}:=(-\varepsilon,\varepsilon)\times(-d,d)$.

To proceed with a proof we need the following auxiliary result [9]:

######  Lemma 4.4.

Let $\Pi^{\prime}\subset\mathbb{R}^{n}$ be a convex set and let $G$ and $Q$ be
arbitrary measurable sets in $\Pi^{\prime}$ with $\mu\,(G)\neq 0$. Then, for
all $v\in\mathcal{H}^{1}(\Pi^{\prime})$, the following inequality holds:

$\displaystyle\int_{Q}|v|^{2}\,dx\,dy$ (4.31)
$\displaystyle\leq\frac{2\mu\,(Q)}{\mu\,(G)}\int_{G}|v|^{2}\,dx\,dy+\frac{C(n)(d(\Pi^{\prime}))^{n+1}(\mu\,(Q))^{1/n}}{\mu\,(G)}\int_{\Pi^{\prime}}|\nabla
v|^{2}\,dx\,dy,$

where $d(\Pi^{\prime})$ is the parameter of $\Pi^{\prime}$, $\mu$ is the
Lebesque measure on $\mathbb{R}^{n}$, and the constant $C(n)$ depends only on
the dimension of $\mathbb{R}^{n}$.

Let $G=\Pi^{\prime}=\Pi(d)$ and $Q=\Omega_{\varepsilon}$. Then using the Lemma
4.4 for the function $f$ we get

$\displaystyle\int_{\Omega_{\varepsilon}}|f|^{2}\,dx\,dy$
$\displaystyle\leq\frac{8\varepsilon
d}{\mu(\Pi(d))}\int_{\Pi(d)}|f|^{2}\,dx\,dy+\frac{2C(2)(d(\Omega))^{3}\sqrt{\varepsilon
d}}{\mu(\Pi(d))}\int_{\Pi(d)}|\nabla f|^{2}\,dx\,dy,$
$\displaystyle\int_{\Omega_{\varepsilon}}|\nabla f|^{2}\,dx\,dy$
$\displaystyle\leq\frac{8\varepsilon d}{\mu(\Pi(d))}\int_{\Pi(d)}|\nabla
f|^{2}\,dx\,dy+\frac{2C(2)(d(\Omega))^{3}\sqrt{\varepsilon
d}}{\mu(\Pi(d))}\int_{\Pi(d)}\left(\left|\nabla\left(\frac{\partial
f}{\partial x}\right)\right|^{2}+\left|\nabla\left(\frac{\partial f}{\partial
y}\right)\right|^{2}\right)\,dx\,dy$ $\displaystyle\leq\frac{8\varepsilon
d}{\mu(\Pi(d))}\int_{\Pi(d)}|\nabla
f|^{2}\,dx\,dy+\frac{2C(2)(d(\Omega))^{3}\sqrt{\varepsilon
d}}{\mu(\Pi(d))}\int_{\Pi(d)}\left(\left|\frac{\partial^{2}f}{\partial
x^{2}}\right|^{2}+\left|\frac{\partial^{2}f}{\partial
y^{2}}\right|^{2}+2\left|\frac{\partial^{2}f}{\partial x\partial
y}\right|^{2}\right)\,dx\,dy\,.$

With the above bounds, one can show that sufficiently small values of
$\varepsilon$, the following is true

$\int_{\Omega_{\varepsilon}}(|\nabla
f|^{2}+|f|^{2})\,dxdy\leq\tilde{c}\,\sqrt{\varepsilon}\,\|f\|_{\mathcal{H}^{2}(\Pi(d))},$
(4.32)

where

$\tilde{c}=\frac{2C(2)(d(\Omega))^{3}\sqrt{d}}{\mu(\Pi(d))}=\frac{C(2)(d(\Omega))^{3}}{2d\sqrt{d}}$

and $\|f\|_{\mathcal{H}^{2}(\Pi(d))}$ means the Sobolev
$\mathcal{H}^{2}(\Pi(d))$ norm of $f$.

Next we need the following interior regularity theorem [1]:

######  Theorem 4.5.

(Interior Regularity Theorem.) Suppose that $h\in\mathcal{H}^{1}(\Omega)$ is a
weak solution of $-\Delta h=w$. Then
$h\in\mathcal{H}^{2}_{\mathrm{loc}}(\Omega)$ and for each
$\Omega_{0}\subset\Omega$ there exists a constant $c=c(\Omega_{0})$
independent of $h$ and $w$ such that:

$\|h\|_{\mathcal{H}^{2}(\Omega_{0})}\leq
c\left(\|h\|_{L^{2}(\Omega)}+\|w\|_{L^{2}(\Omega)}\right).$ (4.33)

In view of the above theorem the right-hand side of (4.32) can be estimated as
follows

$\int_{\Omega_{\varepsilon}}(|\nabla f|^{2}+|f|^{2})\,dx\,dy\leq
c\,\tilde{c}\,\sqrt{\varepsilon}\,\int_{\Omega}(|f|^{2}+|\Delta
f|^{2})\,dx\,dy,$ (4.34)

with some constant $c=c(d)$ does not depend on $\varepsilon$.

By virtue of inequalities (4.30), (4.34) and Lemma 6.3 (see Appendix) we have

$\displaystyle\int_{\partial K_{\varepsilon}}|f(x,y)|^{2}\,d\mu\leq
2c\,\tilde{c}\,\mathrm{max}\left\\{1/d,2d\right\\}\left(1+\|y_{1}^{\prime}\|_{L^{\infty}(-1,1)}^{2}\right)^{1/2}\,\sqrt{\varepsilon}\,\,\int_{\Omega}(|f|^{2}+|\Delta
f|^{2})\,dxdy$ $\displaystyle\leq
2c\,\tilde{c}\,\mathrm{max}\left\\{1/d,2d\right\\}\left(1+\|y_{1}^{\prime}\|_{L^{\infty}(-1,1)}^{2}\right)^{1/2}\,\sqrt{\varepsilon}\,\|f\|_{2}^{2}\,.$

The above combined the fact that

$\gamma_{K_{\varepsilon}}\int_{\partial
K_{\varepsilon}}|u|^{2}\,d\mu\leq\|u\|_{1}^{2}$

estimates the right-hand side of the inequality (4.28) as follows

$\displaystyle\left|\int_{\partial
K_{\varepsilon}}(J_{1}f)\,\overline{u}\,d\mu\right|$
$\displaystyle\leq\left(\frac{2c\,\tilde{c}\,\mathrm{max}\left\\{1/d,2d\right\\}}{\gamma_{K_{\varepsilon}}}\right)^{1/2}\left(1+\|y_{1}^{\prime}\|_{L^{\infty}(-1,1)}^{2}\right)^{1/4}\,\varepsilon^{1/4}\,\|f\|_{2}\|u\|_{1}\,.$

By virtue of (4.26), (4.27) and the above inequality the right-hand side of
inequality (4.20) satisfies

$\mathrm{r.h.s.}(\ref{7''})=\mathcal{O}\left(\varepsilon^{1/6}+\gamma_{K_{\varepsilon}}^{1/2}\varepsilon^{1/4}\right)\,\|f\|_{2}\,\|u\|_{1},\quad\varepsilon\to
0,$

which ends the proof.∎

##  5 Proof of Lemma 4.3

Let $M_{\epsilon}$ be the subset of $\mathbb{B}_{\epsilon}\setminus\Gamma$
such that for every point $(x_{0},y_{0})\in M_{\epsilon}$ the line
$l_{x_{0}}:=\\{x=x_{0}\\}$ intersects the boundary of $\Gamma$ at least once.
Let $(x_{0},y_{1}(x_{0}))$ be the point of intersection of $l_{x_{0}}$ with
$\partial\Gamma$ and let $(x_{0},y_{2}(x_{0}))\in\Omega$ be the point such
that the open interval connecting $(x_{0},y_{1}(x_{0},y_{0}))$ and
$(x_{0},y_{2}(x_{0},y_{0}))$ belongs to $\Omega$ and has an empty intersection
with $\Gamma$. Without loss of generality suppose that
$y_{1}(x_{0})<y_{2}(x_{0})$. Then for any
$u\in\mathcal{H}^{1}(\mathbb{B}_{\varepsilon}\setminus\Gamma)$ and almost all
$(x_{0},y_{0})\in M_{\epsilon}$ we have

$u(x_{0},y_{0})=u(x_{0},y_{1}(x_{0},y_{0}))+\int_{y_{1}(x_{0},y_{0})}^{y}\frac{\partial
u}{\partial t}(x_{0},t)\,dt.$

Let $P_{\epsilon}$ denote the projection of $M_{\epsilon}$ onto the axis
$\mathbb{X}$. We get

$\displaystyle\int_{M_{\epsilon}}|u(x,y)|^{2}\,dx\,dy$ $\displaystyle\leq
2\int_{M_{\epsilon}}|u(x,y_{1}(x,y))|^{2}\,dx\,dy+2\int_{M_{\epsilon}}\left|\int_{y_{1}(x,y)}^{y}\frac{\partial
u}{\partial t}(x,t)\,dt\right|^{2}\,dx\,dy$
$\displaystyle=2\int_{P_{\epsilon}}\int_{\\{y:\,(x,y)\in
M_{\epsilon}\\}}|u(x,y_{1}(x))|^{2}\,dxdy$
$\displaystyle+2\int_{P_{\epsilon}}\int_{\\{y:\,(x,y)\in
M_{\epsilon}\\}}\left|\int_{y_{1}(x,y)}^{y}\frac{\partial u}{\partial
t}(x,t)\,dt\right|^{2}\,dx\,dy$ $\displaystyle\leq
2\int_{P_{\epsilon}}\int_{\\{y:\,(x,y)\in
M_{\epsilon}\\}}|u(x,y_{1}(x))|^{2}\,dx\,dy$
$\displaystyle+2\int_{P_{\epsilon}}\int_{\\{y:\,(x,y)\in
M_{\epsilon}\\}}(y-y_{1}(x))\int_{y_{1}(x)}^{y}\left|\frac{\partial
u}{\partial t}(x,t)\right|^{2}\,dt\,dx\,dy$ $\displaystyle\leq
4\epsilon\int_{P_{\epsilon}}|u(x,y_{1}(x))|^{2}\,dx+8\epsilon^{2}\int_{P_{\epsilon}}\int_{y_{1}(x)}^{y_{2}(x)}\left|\frac{\partial
u}{\partial t}(x,t)\right|^{2}\,dx\,dt$ $\displaystyle\leq
4\epsilon\int_{P_{\epsilon}}|u(x,y_{1}(x))|^{2}\,dx+8\epsilon^{2}\int_{\mathbb{B}_{\epsilon}\setminus\Gamma}\left|\frac{\partial
u}{\partial t}(x,t)\right|^{2}\,dx\,dt$ $\displaystyle\leq
4\epsilon\int_{P_{\epsilon}}|u(x,y_{1}(x))|^{2}\,dx+8\epsilon^{2}\|u\|_{1}^{2}.$
(5.1)

Since the set $\\{(x,y_{1}(x))\\},\,x\in P_{\epsilon}$, is the part of
$\Gamma$ then

$\displaystyle\int_{P_{\epsilon}}|u(x,y_{1}(x))|^{2}\,dx\leq\int_{P_{\epsilon}}\left(1+y_{1}^{\prime}(x)^{2}\right)^{1/2}|u(x,y_{1}(x))|^{2}\,dx$
$\displaystyle=\int_{\Gamma\cap\\{(x,y):x\in
P_{\epsilon}\\}}|u|^{2}\,d\mu\leq\int_{\Gamma}|u|^{2}\,d\mu\,.$

Returning to the inequality (5.1) and combining the above bound together with
the fact that

$\int_{\Gamma}|u|^{2}\,d\mu\leq\frac{1}{\gamma_{\Gamma}}\|u\|_{1}^{2}$

we get

$\int_{M_{\epsilon}}|u(x,y)|^{2}\,dxdy\leq
4\left(\frac{\epsilon}{\gamma_{\Gamma}}+2\epsilon^{2}\right)\|u\|_{1}^{2}\,.$
(5.2)

Now let us go to the subset
$\left(\mathbb{B}_{\epsilon}\setminus\Gamma\right)\setminus M_{\epsilon}$. For
any
$(x_{0},y_{0})\in\left(\mathbb{B}_{\epsilon}\setminus\Gamma\right)\setminus
M_{\epsilon}$ let $(x_{0},y_{3}(x_{0}))$ with $y_{3}(x_{0})<y_{0}$, be a point
of intersection of line $l_{x_{0}}$ with the boundary of $\Omega$. One can
easily check that there is $y_{4}(x_{0})\in(y_{3}(x_{0}),y_{0})$ such that

$|u(x_{0},y_{4})|\leq\frac{1}{\sqrt{y_{0}-y_{3}(x_{0})}}\sqrt{\int_{y_{3}(x_{0})}^{y_{0}}|u(x_{0},z)|^{2}\,dz}.$

Therefore

$\displaystyle|u(x_{0},y_{0})|^{2}=\left|u(x_{0},y_{4}(x_{0}))+\int_{y_{4}(x_{0})}^{y_{0}}\frac{\partial
u}{\partial z}(x_{0},z)\,dz\right|^{2}$ $\displaystyle\leq
2|u(x_{0},y_{4}(x_{0}))|^{2}+2(y_{0}-y_{4}(x_{0}))\int_{y_{4}(x)}^{y}\left|\frac{\partial
u}{\partial z}(x,z)\right|^{2}\,dy$
$\displaystyle\leq\frac{2}{y_{0}-y_{3}(x_{0})}\int_{y_{3}(x_{0})}^{y_{0}}|u(x_{0},z)|^{2}\,dz+2(y_{0}-y_{3}(x_{0}))\int_{y_{3}(x_{0})}^{y_{0}}\left|\frac{\partial
u}{\partial z}(x_{0},z)\right|^{2}\,dz$
$\displaystyle\leq\frac{2}{\mathrm{dist}((x_{0},y_{0}),\partial\Omega)}\int_{y_{3}(x_{0})}^{y_{0}}|u(x_{0},z)|^{2}\,dz+2\mathrm{diam}(\Omega)\int_{y_{3}(x_{0})}^{y_{0}}\left|\frac{\partial
u}{\partial z}(x_{0},z)\right|^{2}\,dz$
$\displaystyle\leq\widetilde{C}\left(\int_{y_{3}(x_{0})}^{y_{0}}|u(x_{0},z)|^{2}\,dz+\int_{y_{3}(x_{0})}^{y_{0}}\left|\frac{\partial
u}{\partial z}(x_{0},z)\right|^{2}\,dz\right)\,,$

where $\mathrm{diam}(\Omega)$ is the diameter of $\Omega$ and
$\mathrm{dist}((x_{0},y_{0}),\partial\Omega)$ is the distance between
$(x_{0},y_{0})$ and the boundary of $\Omega$ and

$\widetilde{C}=2\mathrm{max}\left\\{\frac{1}{\mathrm{dist}((x_{0},y_{0}),\partial\Omega)},\,\mathrm{diam}(\Omega)\right\\}.$

Let $P_{\epsilon}^{\prime}$ be the projection of
$\left(\mathbb{B}_{\epsilon}\backslash\Gamma\right)\setminus M_{\epsilon}$
onto axis $\mathbb{X}$. So we get

$\displaystyle\int_{\left(\mathbb{B}_{\epsilon}\backslash\Gamma\right)\setminus
M_{\epsilon}}|u(x,y)|^{2}\,dxdy\leq\widetilde{C}\int_{\left(\mathbb{B}_{\epsilon}\backslash\Gamma\right)\backslash
M_{\epsilon}}\int_{y_{3}(x)}^{y}\left(\left|\frac{\partial u}{\partial
z}(x,z)\right|^{2}+|u(x,z)|^{2}\right)\,dzdxdy$
$\displaystyle=\widetilde{C}\int_{P^{\prime}_{\epsilon}}\int_{\\{y:(x,y)\in\left(\mathbb{B}_{\epsilon}\backslash\Gamma\right)\backslash
M_{\epsilon}\\}}\int_{y_{3}(x)}^{y}\left(\left|\frac{\partial u}{\partial
z}(x,z)\right|^{2}+|u(x,z)|^{2}\right)\,dzdxdy\,.$

Then finally

$\int_{\left(\mathbb{B}_{\epsilon}\backslash\Gamma\right)\setminus
M_{\epsilon}}|u(x,y)|^{2}\,dxdy\leq
2\widetilde{C}\epsilon\int_{\Omega\backslash\Gamma}(|\nabla
u|^{2}+|u|^{2})\,dxdy\leq 2\widetilde{C}\epsilon\|u\|_{1}^{2}\,.$

This together (5.2) proves the lemma with

$C_{1}=\mathrm{max}\\{2\widetilde{C}+8\varepsilon,4\\}<\mathrm{max}\\{2\widetilde{C}+8,4\\}.$

∎

##  6 Appendix

In this section we mention several useful lemmas proved in [2].

######  Lemma 6.1.

Let $\Omega^{\prime}$ be an open bounded domain in $\mathbb{R}^{2}$ and let
$p\in\Omega^{\prime}$ be some fixed point. Suppose that
$\mathbb{B}_{\epsilon}\subset\Omega^{\prime}$ is a ball with center at $p$ and
radius $\epsilon>0$. For any function $g\in\mathcal{H}^{2}(\Omega^{\prime})$
the following estimate takes place

$\int_{\mathbb{B}_{\epsilon}}|\nabla g|^{2}\,dxdy\leq
C_{2}\epsilon^{4/3}\int_{\Omega^{\prime}}|-\Delta g+g|^{2}\,dxdy,$

with the constant $C_{2}$ depending on $\Omega$.

######  Lemma 6.2.

Let $\mathbb{B}_{2\varepsilon}$ and $\mathbb{B}_{\varepsilon}$ be the balls
centered on the origin and the radii $\varepsilon$ and $2\varepsilon$. Let
$g\in\mathcal{H}^{1}(\mathbb{B}_{2\varepsilon}\setminus\mathbb{B}_{\varepsilon})$.
Then there exists $\tau\in(\varepsilon,2\varepsilon)$ such that

$\int_{0}^{2\pi}\left|\frac{\partial\tilde{g}}{\partial\varphi}(\tau,\varphi)\right|^{2}\,d\varphi\leq
4\int_{\mathbb{B}_{2\varepsilon}\setminus\mathbb{B}_{\varepsilon}}(|\nabla
g|^{2}+|g|^{2})\,dxdy,$

where $\tilde{g}(r,\varphi):=g(r\cos\varphi,r\sin\varphi)$.

######  Lemma 6.3.

For any $z\in\mathrm{Dom}(H_{\Omega}^{N})$ the following estimate is valid

$\int_{\Omega}|-\Delta z+z|^{2}\,dxdy\geq\int_{\Omega}(|\Delta
z|^{2}+|z|^{2})\,dxdy.$

### Acknowledgements

D.B. acknowledges support from the Czech Science Foundation (GACR), the
project 21-07129S.

## References

  * [1] L. Beck, Elliptic Regularity Theory, Springer International Publishing, Switzerland, 2016
  * [2] D. Barseghyan, B. Schneider and L.H. Hai, Neumann Laplacian in a perturbed domain, The Mediterranean Journal of Mathematics (2022), 19(126), 1–17.
  * [3] D.I. Borisov, A.I. Mukhametrakhimova, On norm resolvent convergence for elliptic operators in multi-dimensional domains with small holes, Journal of Mathematical Sciences 232 (2018), 283–298.
  * [4] M. Balzano and L. Notarantonio, On the asymptotic behavior of Dirichlet problems in a Riemannian manifold less small random holes, Rendiconti del Seminario Matematico della Universita di Padova 100 (1998), 249–282.
  * [5] A. Colette, O. Post, Wildly perturbed manifolds: norm resolvent and spectral convergence, Journal of Spectral Theory, European Mathematical Society, In press.
  * [6] D. Daners, Dirichlet problems on varying domains, Journal of Differential Equations 188(2003), 591–624.
  * [7] L. C. Evans, Partial differential equations, American Mathematical Society, 2010.
  * [8] P. Exner, D. Borisov, Approximation of point interactions by geometric perturbations in two-dimensional domains, Bulletin of Mathematical Sciences(2022) 2250003.
  * [9] V. A. Marchenko, E. Ya. Khruslov, Homogenization of partial differential equations, Progress in Mathematical Physics, 46, Birkhäuser, Boston, 2006.
  * [10] V. G. Maz’ya, S. A. Nazarov, B. A. Plamenewskii, Asymptotic Theory of Elliptic Boundary Value Problems in Singularly Perturbed Domains, I, II. Operator Theory: Advances and Applications, vols. 111, 112. Birkhäuser,Basel (2000) (translation of the original in German published by Akademie Verlag 1991)
  * [11] S. Ozawa, Spectra of domains with small spherical Neumann boundary, Proceedings of the Japan Academy, Series A, Mathematical Sciences 58(5) (1982), 190–192.
  * [12] O. Post, Spectral convergence of quasi-one-dimensional spaces, Annales Henri Poincare 7 (2006), 933– 973.
  * [13] J. Rauch, M. Taylor, Potential and scattering theory on wildly perturbed domains, Journal of Functional Analysis 18(1975), 27–59.
  * [14] B. Schweizer, The low-frequency spectrum of small Helmholtz resonators, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences: 20140339, 2015.


# Deep learning inter-atomic potential for irradiation damage in 3C-SiC

Yong Liu  These authors contributed equally to this work  State Key Laboratory
of Nuclear Physics and Technology, School of Physics, CAPT, HEDPS, and IFSA,
College of Engineering, Peking University, Beijing 100871, P. R. China  Hao
Wang  These authors contributed equally to this work  Institute of Nuclear
Physics and Chemistry, China Academy of Engineering Physics, Mianyang 621999,
P. R. China  Linxin Guo  State Key Laboratory of Nuclear Physics and
Technology, School of Physics, CAPT, HEDPS, and IFSA, College of Engineering,
Peking University, Beijing 100871, P. R. China  Zhanfeng Yan  Institute of
Nuclear Physics and Chemistry, China Academy of Engineering Physics, Mianyang
621999, P. R. China  Jian Zheng  Institute of Nuclear Physics and Chemistry,
China Academy of Engineering Physics, Mianyang 621999, P. R. China  Wei Zhou
Institute of Nuclear Physics and Chemistry, China Academy of Engineering
Physics, Mianyang 621999, P. R. China  Jianming Xue  Corresponding author:
<EMAIL_ADDRESS>State Key Laboratory of Nuclear Physics and Technology,
School of Physics, CAPT, HEDPS, and IFSA, College of Engineering, Peking
University, Beijing 100871, P. R. China

###### Abstract

We developed and validated an accurate inter-atomic potential for molecular
dynamics simulation in cubic silicon carbide (3C-SiC) using a deep learning
framework combined with smooth Ziegler-Biersack-Littmark (ZBL) screened
nuclear repulsion potential interpolation. Comparisons of multiple important
properties were made between the deep-learning potential and existing
analytical potentials which are most commonly used in molecular dynamics
simulations of 3C-SiC. Not only for equilibrium properties but also for
significant properties of radiation damage such as defect formation energies
and threshold displacement energies, our deep-learning potential gave closer
predictions to DFT criterion than analytical potentials. The deep-learning
potential framework solved the long-standing dilemma that traditional
empirical potentials currently applied in 3C-SiC radiation damage simulations
gave large disparities with each other and were inconsistent with ab-initio
calculations. A more realistic depiction of the primary irradiation damage
process in 3C-SiC can be given and the accuracy of classical molecular
dynamics simulation for cubic silicon carbide can be expected to the level of
quantum mechanics.

##  I Introduction

Cubic silicon carbide has been widely used for electronic and nuclear
applications due to its outstanding mechanical properties, high thermal
conductivity, chemical stability, and good radiation response [1, 2]. The
mechanical and electrical properties of 3C-SiC are degraded due to the changes
in microstructure when it subjects to high energy neutron in the nuclear
environment. Understanding the primary irradiation process is of crucial
importance to estimate the usable lifetime of this material.

Ab-initio molecular dynamics (AIMD) with density functional theory (DFT) and
classical molecular dynamics (CMD) are the main tools to simulate the primary
irradiation damage process at the atomic level beyond the limits of
experimental techniques [3]. On the one hand, AIMD is accurate but
computationally cost, which can only involve a few hundred atoms and several
hundred picoseconds long [4]. Many thousands of atoms can be knocked out of
equilibrium position by one energetic ion or neutron generated from a nuclear
reaction. Therefore, AIMD can not cover the atomic scale required to simulate
the primary radiation damage process [5]. On the other hand, CMD is efficient
enough to satisfy the computational demand of primary irradiation dynamics
simulation but the accuracy of simulation results greatly depends upon the
employed inter-atomic potential.

Table 1: The widely used empirical potentials for MD simulations of silicon carbide materials[6]. | Potentials | Applications
---|---|---
| Tersoff[7, 8, 9, 10, 11] | Thermal properties, Mechanical properties,
|  | Electrical properties, Polishing,
|  | Ion implantation, Crystal growth,
|  | Irradiation damage, Amorphization,
|  | Fatigue damage, Shock damage
| Tersoff/ZBL[12] | Ion implantation, Irradiation damage
| GW[13] | Irradiation damage
| GW/ZBL[14] | Crystal growth, Irradiation damage
| Vashishta[15] | Mechanical properties, Electrical properties,
|  | Deposition, Shock damage, Polishing
| MEAM[16, 17] | Crystal growth, Thermal properties,
|  | Irradiation damage
| EDIP[18, 19] | Mechanical properties

The widely used potential functions for CMD simulations of silicon carbide
materials and their applications were summarized in Table 1. Although the
expression forms of different empirical analytic potentials are distinguished,
the processes of their development are basically the same. First, a
mathematical function based on a physical understanding of interatomic
interactions in the material was proposed, with a handful of global fitting
parameters. Then a series of labeled physical properties from experimental or
ab-initio calculations were used to fit these adjustable parameters. Finally,
this fixed expression will be used for predicting the energies and forces of
the new configurations in MD simulations. Although all these potential
functions in Table 1 have taken the three-body effect and bond-angle effect
into account to describe the many-body interaction in material and strong
directionality of the covalent bonds, the true interactions in silicon carbide
are determined by complex many-body interactions. The ability of traditional
analytical force fields to fit the corresponding potential energy surface is
inherently limited by their relatively simple functional forms and a few
adjustable parameters. G. Lucas and L. Pizzagalli pointed out that the use of
available empirical potentials is the largest source of errors to calculate
threshold displacement energies in 3C-SiC and called for the improvement of
existing potentials[20]. G.D.Samolyuk’s study shows that the most commonly
used Tersoff and MEAM potentials for SiC are inconsistent with the ab-initio
calculation of defect energetics. Tersoff potential predicts a very high
interstitial formation energy and high defect migration energy[21]. GW-ZBL
potential gives a more realistic description of defect formation energy but
still overestimates the defect migration energy barrier. [21]. Andrey Sarikov
got divergent simulation results from different potentials (including Tersoff,
Vashishta [22] and an analytical bond order potential [11]) in their study of
partial dislocations and stacking faults in 3C-SiC [23]. The inaccurate
depiction of these key physical quantities makes us lose confidence in the
correctness of MD simulation results of radiation in 3C-SiC. A new potential
that can accurately describe the inter-atomic interactions is urgently needed
to be developed.

Recently, machine learning methods combined with DFT training data to build
potential energy surfaces (PES) have been developed rapidly [24, 25, 26, 27,
28, 29]. Compared with the construction method of traditional empirical
potential, machine-learning potentials have a more powerful fitting ability
due to their unintended preset expressions and abundant adjustable parameters
[30]. Moreover, unlike the empirical potentials, which only fit a subset of
properties, machine-learning potentials can sample configurations to train the
PES as many as needed. Due to more general expressions and more complete
training data, machine learning potentials can give a more accurate prediction
of the PES to capture the underlying physical mechanism. A variety of CMD
simulations with DFT accuracy in different areas have been carried out with
the help of machine-learning potentials [31, 32, 33, 34]. In the field of
radiation damage, a set of machine learning potentials have also been
developed to simulate the irradiation damage processes for different materials
such as fcc-aluminum[35], tungsten [29], silicon[36], and bcc-iron [37]. So
far, most of the machine learning potentials are for single-substance systems
because the number of configurations needed to train the model increases
exponentially with the increase of the number of principle elements.

In this work, we applied the DP-ZBL (Deep-learning Potential interpolated with
ZBL) model [35] to train a deep-learning inter-atomic potential hybrid with
ZBL screened nuclear repulsion potential for 3C-SiC. In order to capture the
right pictures when atoms are extremely close to each other, which are
frequent events happening in the irradiation damage process due to the high
kinetic energy of atoms, the generally used ZBL screened nuclear repulsion
potential [38] has been interpolated into the deep learning framework so that
short-range repulsion interaction between atoms can be accurately described.
Here we refer Ref[35] for more details about the interpolation mechanism.
Compared to the analytics empirical potential including Tersoff, MEAM,
Vashishta, EIDP, and GW-ZBL, the DP-ZBL potential not only get the DFT
accuracy for the near-equilibrium properties such as lattice constant, elastic
properties, equation of state, and phonon dispersion but also give a correct
description of short-range repulsion interaction. We put our concentration on
correct prediction for defect formation energies and threshold displacement
energies because these physical quantities play decisive roles in the
irradiation process and there are large disparities between different studies
using distinct inter-atomic potentials. The DP-ZBL model can terminate these
controversies and get a more realistic molecular dynamics simulation of
radiation damage in 3C-SiC.

##  II Method

###  II.1 Training process

The accuracy and transferability of DP models are determined by the quality of
the training dataset. The training dataset should be complete enough to cover
the target simulation space. To get a good description of the energy and force
of dimer (Si-Si, Si-C, C-C), elastic properties, phonon dispersion, and defect
formation energies, the corresponding configurations recorded in Table 2 are
taken as the initial training dataset.

Data used to train the neural network were all generated by DFT calculations
with VASP code [39, 40, 41]. General gradient approximate (GGA) with PBE [42]
exchange-correlation functional has been used. The plane-wave cutoff energy
was set high enough to 600eV to cover large deformation in the irradiation
process. Consistent spacing between k-points in Brillouin zone (KSPACING =
0.15 $\rm{\AA^{-1}}$) was integrated using Gamma centered grid for all
configurations. Gaussian smearing with 0.03 eV width was applied to help
convergence. Spin-polarized calculations were considered to account for the
possible spin polarization of various defect configurations.

Table 2: Configurations included in the initial training database for corresponding properties. | Concerned properties | Configuration type
---|---|---
| Bulk properties | Equilibrium state
|  | Compressed
|  | Stretched
| Thermo properties | Atom displaced
| Elastic properties | Elastically distorted
| Defect properties | Vacancy with strain
|  | Antisite with strain
|  | Tetrahedral interstitial with strain
| Liquid phase | Frames of liquid trajectory
| Short-range interactions | Dimers
| Irradiation damage | Frames of PKA activation trajectory

After the initial training dataset was used to kick off the first round train,
an active learning training process was performed by DPGEN[43] to sample more
configurations into the training dataset.The active learning process
terminated when we validated that the potential energy is good enough. In this
work, we went through fifty active learning iterations to get the final
potential energy model, which made the training dataset sampling 33898
configurations in total. The temperature of NVT setting in the exploration
stage went up from 300K to 4000K and the environmental pressure is set to one
atmosphere.

###  II.2 Interpolation for short range repulsion

We refer to our previous work Ref[35] to get the details of the interpolation
method between DP and ZBL. Meanwhile, R.E.Stoller’s work in Ref[44] provides a
more systematic and effective procedure to bridge the equilibrium and short-
range parts, which can avoid lots of invalid tests. First, the table recording
the energy of dimer Si-Si, Si-C, C-C at a short range
($\rm{0.001\AA,1.200\AA}$) with $\rm{0.001\AA}$ step was generated by ZBL
formula. The dimer configurations of Si-Si, Si-C, and C-C in range
($\rm{0.5,5.0}$) with $\rm{0.05\AA}$ step were calculated by DFT and added to
the training dataset. Then the DP and the ZBL potential were smoothly docked
in the interval ($\rm{1.0\AA,1.2\AA}$) to let ZBL plays its role in the short
range and the DP model dominates at equilibrium condition.

As shown in Figure 1, when Si-C dimer’s distance is less than 1.0 $\rm{\AA}$,
DP-ZBL is consistent with pure ZBL potential and when Si-C dimer’s distance is
larger than 1.2 $\rm{\AA}$, DP-ZBL gives results dominated by DFT calculation.
And for the Si-C dimer’s distance between 1.0 and 1.2 $\rm{\AA}$, DP-ZBL has
smoothly switched from ZBL to DFT.

Figure 1: The calculated energy versus distance curve of Si-C dimer by DFT,
ZBL and DP-ZBL.

##  III Result

We compared the energies calculated by the DP-ZBL model and the DFT method for
all the configurations in the training dataset. As shown in Fig. 2 and Fig. 3,
DP-ZBL prediction is consistent with DFT calculation as the points basically
distribute around the y=x reference line. The root-mean-squared-errors (RMSEs)
of the energies and the forces are 0.01 $\rm{meV/atom}$ and 0.16 $\rm{eV/\AA}$
respectively, which are within the accuracy allowed for typical Deepmd-kit
training [45] compared with the range of energy and force.

Figure 2: Comparison of energy prediction by DFT and DP for the final training
set. Both axes represent the energy of the configuration divided by the number
of total atoms in the configuration. Figure 3: Comparison of force prediction
by DFT and DP for the final training set. The color bar indicates the density
of data. Force greater than 50 eV/$\rm{\AA}$ are not shown in the figure.

Several near-equilibrium material properties have been calculated by the DP-
ZBL and the DFT method and four different empirical potentials including
Tersoff-ZBL, GW-ZBL, MEAM, and Vashishta for comparison, as shown in Table 3.
including lattice parameter $\rm{a_{0}}$, elastic constants $\rm{C_{11}}$,
$\rm{C_{12}}$, $\rm{C_{44}}$, bulk modulus $\rm{K_{H}}$, Young’s modulus
$\rm{E_{H}}$, shear modulus $\rm{G_{H}}$ (all the moduli are in Hill noting),
cohesive energy $\rm{E_{coh}}$, and defect formation energy. All results were
in excellent agreement with DFT data computed in this work or referred from
other works. In contrast, most empirical potentials do not give good
predictions for most properties. For instance, the Tersoff-ZBL potential
underestimates the lattice constant which is the most basic physical quantity
in simulations. However, Tersoff-ZBL works well with the elastic contants
including $\rm{C_{11}}$, $\rm{C_{12}}$, $\rm{C_{44}}$. The GW-ZBL potential
almost mispredicts the elastic response relationship. Meanwhile, the MEAM
potential has a slight error and the Vashishta potential underestimates the
$\rm{C_{44}}$ constant.

The formation energy of defects in 3C-SiC is defined as follows:

$\rm{E_{f}=E_{defect}-E_{perfect}+n_{Si}\mu_{Si}+n_{C}\mu_{C}}$ (1)

where $\rm{E_{perfect}}$ and $\rm{E_{defect}}$ are the total energy of a
perfect 3C-SiC supercell and a supercell containing a defect, respectively.
The integer $\rm{n_{Si}}$ gives the number of Si atoms removed from
($\rm{n_{Si}\textgreater 0}$) or add to ($\rm{n_{Si}\textless 0}$) the perfect
supercell, and $\rm{n_{C}}$ follows the same logic. The $\rm{\mu_{Si}}$ and
$\rm{\mu_{C}}$ are respectively the chemical potential of the Si atom and the
C atom in the 3C-SiC bulk environment. In this work, all defect formation
energies were calculated for the Si-rich condition, which means the chemical
potential of the Si atom in 3C-SiC is limited to that in the cubic silicon
crystal. In this context, $\rm{\mu_{Si}=\mu_{Si}(bulk)}$ and
$\rm{\mu_{C}=\mu_{SiC}-\mu_{Si}}$, where $\rm{\mu_{SiC}}$ is the chemical
potential of Si-C atom pair in 3C-SiC [46]. The results calculated by the DP-
ZBL model match well with the DFT results for all the defect configurations in
this work. The GW-ZBL potential underestimates the defect formation energies
of most configurations except the antisite of C replacing Si. The Tersoff-ZBL
potential underestimates the defect formation energies of vacancies and
antisites but overestimates the defect formation energies of the tetrahedral
interstitial. Defect formation energy is a key physical quantity that reflects
the accuracy of irradiation simulation. To sum up, the DP-ZBL model performs
much better than the four listed empirical potentials for the predictions of
the near-equilibrium properties.

Table 3:  Basic properties of 3C-SiC: lattice constant $\rm{a}$, cohesive energy $\rm{E_{coh}}$, elastic constants $C_{\rm{11}}$, $C_{\rm{12}}$ and $C_{\rm{44}}$, bulk modulus $K_{\rm{H}}$ (Hill), Young’s modulus $E_{\rm{H}}$ (Hill), shear modulus $G_{\rm{H}}$ (Hill), Poisson’s ratio $\rm{v_{H}}$(Hill). | Properties | $\rm{DFT_{Ref}}$ | $\rm{DFT_{Current}}$ | DP-ZBL | Tersoff-ZBL | GW-ZBL | MEAM | Vashishta | EDIP
---|---|---|---|---|---|---|---|---|---
| $\rm{a_{0}(\AA)}$ | 4.3805111Reference [47] | 4.3784 | 4.3778 | 4.2796 | 4.3600 | 4.3595 | 4.3582 | 4.3624
| $\rm{E_{coh}}$ | -15.062411footnotemark: 1 | -15.0643 | -15.0630 | -12.68 | -12.82 | -12.86 | -12.68 | -12.67
| $\rm{C_{11}(GPa)}$ | 383.911footnotemark: 1 | 380.8 | 375.3 | 445.7 | 265.2 | 396.5 | 390.1 | 396.8
| $\rm{C_{12}(GPa)}$ | 127.611footnotemark: 1 | 126.8 | 124.0 | 138.7 | 219.3 | 147.1 | 142.8 | 140.5
| $\rm{C_{44}(GPa)}$ | 239.511footnotemark: 1 | 240.1 | 223.1 | 220.0 | 101.1 | 135.6 | 136.9 | 170.3
| $\rm{K_{H}(GPa)}$ | 213.011footnotemark: 1 | 211.5 | 207.7 | 241.0 | 234.6 | 230.3 | 225.2 | 226.0
| $\rm{E_{H}(GPa)}$ | 432.811footnotemark: 1 | 431.4 | 413.9 | 452.2 | 156.5 | 330.6 | 330.1 | 372.5
| $\rm{G_{H}(GPa)}$ | 186.311footnotemark: 1 | 186.0 | 177.2 | 190.4 | 56.4 | 131.1 | 131.4 | 152.0
| $\rm{v_{H}}$ | 0.1611footnotemark: 1 | 0.16 | 0.17 | 0.19 | 0.39 | 0.26 | 0.26 | 0.23
| $\rm{V_{Si}}$ | 7.75222Reference [48] | 7.72 | 7.66 | 8.24 | 6.89 | 4.90 | 12.73 | 4.60
| $\rm{V_{C}}$ | 4.0922footnotemark: 2 | 4.21 | 4.10 | 3.76 | -0.84 | 1.06 | -3.38 | 1.22
| $\rm{C_{Si}}$ | 3.9422footnotemark: 2 | 3.92 | 3.91 | 3.29 | 8.87 | 2.74 | 33.48 | 3.02
| $\rm{Si_{C}}$ | 3.2922footnotemark: 2 | 3.35 | 3.31 | 4.90 | 0.74 | 3.84 | -3.32 | 2.04
| $\rm{Si_{TSi}}$ | 10.87333Reference [49] | 10.22 | 10.26 | 16.65 | 3.23 | 4.00 | -2.07 | 11.69
| $\rm{Si_{TC}}$ | 9.0422footnotemark: 2 | 8.47 | 8.34 | 16.48 | 0.33 | 3.22 | -3.41 | 12.24
| $\rm{C_{TSi}}$ | 10.0933footnotemark: 3 | 9.97 | 9.88 | 4.89 | 7.86 | 9.08 | 17.84 | 6.69
| $\rm{C_{TC}}$ | 11.1022footnotemark: 2 | 10.96 | 10.86 | 7.89 | 8.22 | 3.05 | 21.16 | 8.29

The equation of state curves of the 3C-SiC phase computed by different
empirical potentials, the DP-ZBL model and the DFT method are illustrated in
Figure 4. The DP-ZBL model well reproduces the DFT results, which indicates
that the DP-ZBL potential is capable to cover the high compressing and
stretching conditions. By contrast, the Tersoff-ZBL and Vashishta potentials
produce large errors compared with the DFT criterion in the compressing
condition. The GW-ZBL and Tersoff-ZBL potentials overestimate the potential
energy when the system is stretched to 1.2 times relative to the equilibrium
state.

Figure 4: Equation of state of the 3C-SiC phase as computed by different
empirical potentials and the DP-ZBL model and the DFT method.

The phonon dispersion curve of the 3C-SiC phase has been calculated along the
high symmetry directions of $\rm{\Gamma}$-X-K-$\rm{\Gamma}$-L by our DP-ZBL
model and the DFT method. The force constants were calculated by density
functional perturbation theory using VASP for the DFT method and were
calculated using PhonoLammps for the DP-ZBL model. Then the Phonopy package
[50] was used to compute the phonon dispersion relations. Non-metallic
crystals are polarized by atomic displacements and the generated macroscopic
field changes force constants near $\rm{\Gamma}$ point[51]. To take this into
consideration, phonon frequencies at general q-points with long-range dipole-
dipole interaction were calculated by the method of Gonze et al. [52, 53]. The
Bron effective charges ($\rm{Z^{*})}$ and dieletric constant ($\epsilon_{0}$)
calculated by GGA functional ($\rm{Z^{*}_{Si}}=\rm{Z^{*}_{C}}=2.69$,
$\epsilon_{0}=6.99$) are in agreement with the experimental value
($\rm{Z^{*}_{Si}}=\rm{Z^{*}_{C}}=2.69$ [54], $\epsilon_{0}=6.52$ [55]). As
shown in Fig. 5, both acoustic branches and optical branches of the six phonon
modes generated by the DP-ZBL model are close to the DFT results. In addition,
the results from both theoretical models match well with the experimental data
measured by Serrano et al. [56] at room temperature using inelastic x-ray
scattering (IXS) with the synchrotron radiation source, which means a good
description of crystal thermal response of 3C-SiC can be predicted with our
DP-ZBL potential.

Figure 5: Phonon dispersion curve of 3C-SiC calculated by our DP-ZBL model and
the DFT method, as well as the experimental data measured by Serrano et al.
[56] using IXS.

Threshold displacement energy ($\rm{E_{d}}$) is defined as the minimum kinetic
energy transferred to a lattice atom to displace it from its original Wigner-
Seitz cell and form a stable Frenkel pair [57]. $\rm{E_{d}}$ is a critical
physical parameter for estimating damage production and predicting the defect
profile under ion, neutron, and electron irradiation [58]. For example,
$\rm{E_{d}}$ is a key input in large-scale irradiation simulation packages
such as SRIM and TRIM to determine implantation profiles in doping processes
or calculate damage accumulation in materials [59]. In this work, the
$\rm{E_{d}}$ for both Si and C primary knock-on atoms (PKAs) along four
typical low-index crystallographic directions including [100], [110], [111],
and $\rm{[\overline{111}]}$ were calculated using different interatomic
potentials for comparison. The simulations were performed at 300 K. A noncubic
simulation box of $\rm{10\times 10\times 12}$ supercell (9600 atoms) with
periodic boundary conditions was used. Kinetic energies in 0.5 eV increments
were progressively assigned to a specific PKA atom in the central area to find
the minimum energy. The simulation system was relaxed in the canonical
ensemble (NVT) for 10 ps at 300K followed by cascade simulation in the
microcanonical ensemble (NVE) for 10 ps. The Wigner-Seitz defect method was
used to identify defects. From the calculation results summarized in Table 4,
the ${\rm{E_{d}}}$ values generated by our DP-ZBL are close to the DFT
calculations carried out by Zhao et al[60]. GW-ZBL, MEAM, and Vashishta show
divergence from the DFT calculation in multiple crystal directions. Tersoff
performs better than the three other empirical potentials, but it is also out
of line with the DFT values for Si PKA in direction $\rm{[\overline{111}]}$
and C PKA in direction [110]. After the comparison, it is clear that the DP-
ZBL potential yield better $\rm{E_{d}}$ values than the empirical potential
functions.

Table 4:  Threshold displacement energy calculated by DFT, DP-ZBL, and a range of empirical interatomic potentials. |  | DFT[60] | DP-ZBL | Tersoff-ZBL | GW-ZBL | MEAM | Vashishata | EDIP
---|---|---|---|---|---|---|---|---
| Si $\rm{[100]}$ | 41 | 33.5 | 47.0 | 20.5 | 36.5 | 29.5 | 42.0
| Si $\rm{[110]}$ | 50 | 47.0 | 41.0 | 23.5 | 26.5 | 22.5 | 42.0
| Si $\rm{[111]}$ | 21 | 23.0 | 26.0 | 12.5 | 25.5 | 46.0 | 21.5
| Si $\rm{[\overline{111}]}$ | 33 | 43.5 | 44.0 | 31.5 | 26.0 | 35.0 | 22.5
| C $\rm{[100]}$ | 18 | 15.0 | 15.5 | 11.5 | 23.0 | 47.5 | 15.5
| C $\rm{[110]}$ | 19 | 17.0 | 29.0 | 15.5 | 27.0 | 47.5 | 15.5
| C $\rm{[111]}$ | 17 | 15.5 | 23.5 | 8.0 | 19.5 | 150.5 | 15.0
| C $\rm{[\overline{111}]}$ | 50 | 43.5 | 53.0 | 23.0 | 33.5 | 15.5 | 16.0

Then we carried out cascade simulations involving $60\times 60\times 60$ unit
cells (1728000 atoms) using different potentials including Tersoff-ZBL, MEAM,
GW-ZBL, and our DP-ZBL, which are widely used in the simulations of
irradiation damage. The simulated system was equilibrated for 10 ps with
timesteps of 1 fs in NVT ensemble at 300 K. Then a single Si atom in the
central area was given kinetic energy of 5.0 keV in $\rm{[135]}$ direction to
initialize the cascade while holding zero total momentum of the system. The
cascade evolved for 10 ps in the NVE ensemble and during this period the
timesteps were modified in order that the distance traveled by the fastest
particle in the system was less than 0.1 $\rm{\AA}$ per timestep. To dissipate
the heat generated by the PKA, the NVT ensemble of 300K was applied to the
boundary region (2 times the lattice constant, about 8.8 $\rm{\AA}$). The
Wigner-Seitz cell analysis method was used to determine the defects number
with Ovito software [61]. The time-dependent curves of defect amount are shown
in Figure 6. The peak of defects amount with the GW-ZBL potential is much
higher than the results calculated by the other three potentials and so is the
number of residual defects. We infer this is because the GW-ZBL potential
underestimates the threshold displacement energy of either silicon or carbon
atom. There are more than 65% recombine of interstitials and vacancies during
the annealing process with the GW-ZBL potential and the DP-ZBL potential,
while only 20% $\rm{\sim}$ 30% recombine with the two other potentials. Defect
yields and ratios predicted by the four potentials are distinguishing. Among
the different point defects, vacancies and interstitials of carbon and silicon
are dominant whether in thermal peak or stable state as shown in Figure 7. The
DP-ZBL model predicts slightly more vacancies and interstitials of silicon
than carbon at the thermal peak and residual vacancies and interstitials of
carbon but domains after annealing.

Figure 6: Comparison of the production of vacancies time dependent
relationship caused by a single 5 keV Si PKA. Figure 7: Peak (slightly
transparent) and stable numbers (solid) of defects caused by a single 5 keV Si
PKA.

##  IV Conclusions

In this work, a potential energy surface for silicon carbide was developed
with our DP-ZBL model. Compared with the four most commonly used empirical
interatomic potentials for SiC, the DP-ZBL potential can not only give a
better performance on the prediction of near-equilibrium properties including
lattice constant, elastic coefficients, equation of state, phonon dispersion,
and defect formation energies but also depict a more precise picture of
irradiation damage. More accurate values of key parameters in irradiation such
as threshold displacement energy and defect migration energy can be gotten by
using the DP-ZBL potential. Furthermore, our work provides a feasible approach
to figuring out the primary irradiation damage process in covalent compound
materials with ab-initio accuracy.

###### Acknowledgements.

This work is supported by National Natural Science Foundation of China (Grants
No.12135002, No.12205269 and No.U20B2010), fund of Science and Technology on
Plasma Physics Laboratory (No.22ZZJJ0601) and the Nuclear Energy Development
Project. We are grateful for computing resources provided by High performance
Computing Platform of Peking University.

## References

  * Gao  _et al._ [2004] F. Gao, W. J. Weber, M. Posselt, and V. Belko, Physical Review B 69, 245205 (2004).
  * Oda  _et al._ [2013] T. Oda, Y. Zhang, and W. J. Weber, The Journal of chemical physics 139, 124707 (2013).
  * Nordlund [2019] K. Nordlund, Journal of Nuclear Materials 520, 273 (2019), 00066.
  * Zhang  _et al._ [2018a] L. Zhang, J. Han, H. Wang, R. Car, and E. Weinan, Physical review letters 120, 143001 (2018a).
  * Dudarev [2013] S. Dudarev, Annual Review of Materials Research 43, 35 (2013).
  * Yan  _et al._ [2023] Z. Yan, R. Liu, B. Liu, Y. Shao, and M. Liu, Energies 16, 1176 (2023).
  * Tersoff [1989] J. Tersoff, Physical review B 39, 5566 (1989).
  * Halicioglu [1995] T. Halicioglu, Physical Review B 51, 7217 (1995).
  * Tersoff [1994] J. Tersoff, Physical Review B 49, 16349 (1994).
  * Tersoff [1990] J. Tersoff, Physical review letters 64, 1757 (1990).
  * Erhart and Albe [2005] P. Erhart and K. Albe, Physical Review B 71, 035211 (2005).
  * Devanathan  _et al._ [1998] R. Devanathan, T. D. de la Rubia, and W. J. Weber, Journal of nuclear materials 253, 47 (1998).
  * Gao and Weber [2002] F. Gao and W. J. Weber, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 191, 504 (2002).
  * Gao  _et al._ [2001] F. Gao, E. J. Bylaska, W. J. Weber, and L. R. Corrales, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 180, 286 (2001).
  * Vashishta  _et al._ [2007a] P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino, Journal of Applied Physics 101, 103515 (2007a).
  * [16] H. Huangts, N. M. Ghoniemt, J. K. Wongt, and M. I. Baskest, .
  * Kang  _et al._ [2014] K.-H. Kang, T. Eun, M.-C. Jun, and B.-J. Lee, Journal of crystal growth 389, 120 (2014).
  * Bazant  _et al._ [1997] M. Z. Bazant, E. Kaxiras, and J. F. Justo, Physical Review B 56, 8542 (1997).
  * Jiang  _et al._ [2012] C. Jiang, D. Morgan, and I. Szlufarska, Physical Review B 86, 144118 (2012).
  * Lucas and Pizzagalli [2005] G. Lucas and L. Pizzagalli, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 229, 359 (2005).
  * Samolyuk  _et al._ [2015] G. D. Samolyuk, Y. Osetsky, and R. E. Stoller, Journal of Nuclear Materials 465, 83 (2015).
  * Vashishta  _et al._ [2007b] P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino, Journal of applied physics 101, 103515 (2007b).
  * Sarikov  _et al._ [2019] A. Sarikov, A. Marzegalli, L. Barbisan, E. Scalise, F. Montalenti, and L. Miglio, Modelling and Simulation in Materials Science and Engineering 28, 015002 (2019).
  * Behler [2011] J. Behler, The Journal of chemical physics 134, 074106 (2011).
  * Zhang  _et al._ [2018b] Y.-J. Zhang, A. Khorshidi, G. Kastlunger, and A. A. Peterson, The Journal of chemical physics 148, 241740 (2018b).
  * Behler [2016] J. Behler, The Journal of chemical physics 145, 170901 (2016).
  * Botu  _et al._ [2017] V. Botu, R. Batra, J. Chapman, and R. Ramprasad, The Journal of Physical Chemistry C 121, 511 (2017).
  * Zhang  _et al._ [2019] L. Zhang, D.-Y. Lin, H. Wang, R. Car, and E. Weinan, Physical Review Materials 3, 023804 (2019).
  * Byggmästar  _et al._ [2019] J. Byggmästar, A. Hamedani, K. Nordlund, and F. Djurabekova, Physical Review B 100, 144105 (2019).
  * Zhang  _et al._ [2018c] L. Zhang, J. Han, H. Wang, W. A. Saidi, R. Car, _et al._ , arXiv preprint arXiv:1805.09003 (2018c).
  * Niu  _et al._ [2020] H. Niu, L. Bonati, P. M. Piaggi, and M. Parrinello, Nature Communications 11, 2654 (2020), 00000.
  * Zeng  _et al._ [2020a] J. Zeng, L. Cao, M. Xu, T. Zhu, and J. Z. H. Zhang, Nature Communications 11, 5713 (2020a), 00000.
  * Zhang  _et al._ [2021] L. Zhang, H. Wang, R. Car, and W. E, Physical Review Letters 126, 236001 (2021), 00000.
  * Wang  _et al._ [2022a] J. Wang, H. Shen, R. Yang, K. Xie, C. Zhang, L. Chen, K.-M. Ho, C.-Z. Wang, and S. Wang, Carbon 186, 1 (2022a), 00000.
  * Wang  _et al._ [2019] H. Wang, X. Guo, L. Zhang, H. Wang, and J. Xue, Applied Physics Letters 114, 244101 (2019).
  * Hamedani  _et al._ [2021] A. Hamedani, J. Byggmästar, F. Djurabekova, G. Alahyarizadeh, R. Ghaderi, A. Minuchehr, and K. Nordlund, Physical Review Materials 5, 114603 (2021).
  * Wang  _et al._ [2022b] Y. Wang, J. Liu, J. Li, J. Mei, Z. Li, W. Lai, and F. Xue, Computational Materials Science 202, 110960 (2022b).
  * Ziegler and Biersack [1985] J. F. Ziegler and J. P. Biersack, in  _Treatise on heavy-ion science_ (Springer, 1985) pp. 93–129.
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Physical review B 47, 558 (1993).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Physical Review B 49, 14251 (1994).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Computational materials science 6, 15 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Zhang  _et al._ [2020] Y. Zhang, H. Wang, W. Chen, J. Zeng, L. Zhang, H. Wang, and W. E, Computer Physics Communications 253, 107206 (2020).
  * Stoller  _et al._ [2016] R. E. Stoller, A. Tamm, L. K. Béland, G. D. Samolyuk, G. M. Stocks, A. Caro, L. V. Slipchenko, Y. N. Osetsky, A. Aabloo, M. Klintenberg, and Y. Wang, Journal of Chemical Theory and Computation 12, 2871 (2016).
  * Zeng  _et al._ [2020b] J. Zeng, L. Zhang, H. Wang, and T. Zhu, Energy & Fuels 35, 762 (2020b).
  * Chen  _et al._ [2018] N. Chen, Q. Peng, Z. Jiao, I. van Rooyen, W. F. Skerjanc, and F. Gao, Journal of Nuclear Materials 510, 596 (2018).
  * Pizzagalli [2021] L. Pizzagalli, Philosophical Magazine Letters 101, 242 (2021), 00000.
  * Hu  _et al._ [2014] S. Hu, W. Setyawan, R. M. Van Ginhoven, W. Jiang, C. H. Henager, and R. J. Kurtz, Journal of Nuclear Materials 448, 121 (2014), 00008.
  * Sun  _et al._ [2017] D. Sun, R. Li, J. Ding, P. Zhang, Y. Wang, and J. Zhao, Journal of Applied Physics 121, 225111 (2017), 00023.
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, Scripta Materialia 108, 1 (2015).
  * Pick  _et al._ [1970] R. M. Pick, M. H. Cohen, and R. M. Martin, Physical Review B 1, 910 (1970).
  * Gonze  _et al._ [1994] X. Gonze, J.-C. Charlier, D. Allan, and M. Teter, Physical Review B 50, 13035 (1994).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Physical Review B 55, 10355 (1997).
  * Olego  _et al._ [1982] D. Olego, M. Cardona, and P. Vogl, Physical Review B 25, 3878 (1982), 00000.
  * Madelung [1982] O. Madelung, Landolt-Bornstein, Numerical Data and Functional Relationships in Science and Technology, New Series 17 (1982).
  * Serrano  _et al._ [2002] J. Serrano, J. Strempfer, M. Cardona, M. Schwoerer-Böhning, H. Requardt, M. Lorenzen, B. Stojetz, P. Pavone, and W. J. Choyke, Applied Physics Letters 80, 4360 (2002), 00067.
  * Nordlund  _et al._ [2006] K. Nordlund, J. Wallenius, and L. Malerba, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 246, 322 (2006).
  * Wang  _et al._ [2013] X. Wang, H. Y. Xiao, X. Zu, Y. Zhang, and W. J. Weber, Journal of Materials Chemistry C 1, 1665 (2013).
  * [59] G. Lucas and L. Pizzagalli, PHYSICAL REVIEW B , 4.
  * Zhao  _et al._ [2012] S. Zhao, J. Xue, C. Lan, L. Sun, Y. Wang, and S. Yan, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 286, 119 (2012).
  * Stukowski [2009] A. Stukowski, Modelling and simulation in materials science and engineering 18, 015012 (2009).