Proof. Suppose that three functions \(\boldsymbol{\varphi}_1, \boldsymbol{\varphi}_2, \boldsymbol{\varphi}_3\in \ker \mathbf{S}_D\). We first assume that \(\int_{\partial D} \boldsymbol{\varphi}_1 \neq \int_{\partial D} \boldsymbol{\varphi}_2\). Then we can find two constants \(c_1\) and \(c_2\) such that \[\int_{\partial D} \boldsymbol{\varphi}_3 = c_1 \int_{\partial D} \boldsymbol{\varphi}_1 + c_2 \int_{\partial D} \boldsymbol{\varphi}_2.\] We define a new function \[\boldsymbol{\varphi}= \boldsymbol{\varphi}_3 - c_1 \boldsymbol{\varphi}_1 - c_2 \boldsymbol{\varphi}_2.\] It is obvious that \(\mathbf{S}_D[\boldsymbol{\varphi}]=0\) and \(\int_{\partial D} \boldsymbol{\varphi}=0\). From Lemma 2 , we can obtain \(\boldsymbol{\varphi}=0\). If \(\int_{\partial D} \boldsymbol{\varphi}_1 = \int_{\partial D} \boldsymbol{\varphi}_2\), following the same argument, we can obtain that \(\boldsymbol{\varphi}_1 = \boldsymbol{\varphi}_2\). The proof is completed. ◻

Proof of Theorem 1. First, it is noted that the operator \(\hat{\mathbf{S}}_D^{\omega}\) is a Fredholm operator [59]. Thus, we just need to study the kernel of the operator \(\hat{\mathbf{S}}_D^{\omega}\). Suppose that there exists a nontrivial function \(\boldsymbol{\varphi}\in L^2(\partial D)^2\) such that \[\label{eq:kset} \hat{\mathbf{S}}_D^{\omega}[\boldsymbol{\varphi}]=\mathbf{S}_D[\boldsymbol{\varphi}]+\hat{\eta}_{\omega} \int_{\partial D}\boldsymbol{\varphi}=0.\tag{18}\]

\(\bullet\) We first consider the case that the operator \(\mathbf{S}_D\) has a nontrivial kernel. From Lemma 3, we have that the dimension of the kernel of the operator \(\mathbf{S}_D\) is at most \(2\). Thus, we consider the following two cases.

1. If the dimension of the kernel of the operator \(\mathbf{S}_D\) is \(2\), then we can find two functions \(\boldsymbol{\phi}_1, \boldsymbol{\phi}_2 \in \ker \mathbf{S}_D\) such that \(\int_{\partial D} \boldsymbol{\phi}_1 \neq \int_{\partial D} \boldsymbol{\phi}_2\). Multiplying \(\boldsymbol{\phi}_i\), \(i=1,2\), on both sides of 18 and using the fact that the operator \(\mathbf{S}_D\) is self-adjoint, we have that \[\label{eq:kset2} \hat{\eta}_{\omega} \int_{\partial D} \boldsymbol{\varphi}(y) \cdot \int_{\partial D}\boldsymbol{\phi}_i(y) ds(y) =0, \quad i=1,2.\tag{19}\] Then from the choice of \(\omega\) and \(\int_{\partial D} \boldsymbol{\phi}_1 \neq \int_{\partial D} \boldsymbol{\phi}_2\), we obtain \(\int_{\partial D} \boldsymbol{\varphi}=0\). This together with Lemma 2 and the equation 18 yields \(\boldsymbol{\varphi}=0\).

2. If the dimension of the kernel of the operator \(\mathbf{S}_D\) is \(1\), then we can find a function \(\boldsymbol{\phi}\in \ker \mathbf{S}_D\) with \(\int_{\partial D} \boldsymbol{\phi}\neq 0\). Similar to equation 19 , we have that \[\label{eq:kset3} \hat{\eta}_{\omega} \int_{\partial D} \boldsymbol{\varphi}(y) \cdot \int_{\partial D}\boldsymbol{\phi}(y) ds(y) =0.\tag{20}\] Next, we show that the function \(\boldsymbol{\varphi}\) is a real function. In fact, if the function \(\boldsymbol{\varphi}\) is a complex function, then we can write \(\boldsymbol{\varphi}= \boldsymbol{\varphi}_1 + \mathrm{i} \boldsymbol{\varphi}_2\) with \(\boldsymbol{\varphi}_1, \boldsymbol{\varphi}_2\) being real functions. From equation 20 , we obtain that \[\label{eq:kset4} \int_{\partial D} \boldsymbol{\varphi}_2(y) =s_1 \mathbf{a},\tag{21}\] where \(\mathbf{a}= \int_{\partial D} \boldsymbol{\varphi}_1(y)\). Without loss of generality, we assume that \(\int_{\partial D} \boldsymbol{\varphi}_1 \neq 0\) and \(s_1\in \mathbb{R}\). We first consider the case that \(s_1\neq 0\). Substituting the last equation 21 into equation 18 , we have that \[\mathbf{S}_D[\boldsymbol{\varphi}_1] + \mathrm{i} \mathbf{S}_D[\boldsymbol{\varphi}_2] + (\Re\hat{\eta}_{\omega} + \mathrm{i} \Im\hat{\eta}_{\omega}) (1 + \mathrm{i} s_1) \int_{\partial D} \boldsymbol{\varphi}_1 =0.\] Taking the real part and the imaginary part of the last equation, we have that \[\begin{cases} \mathbf{S}_D[\boldsymbol{\varphi}_1] + (\Re\hat{\eta}_{\omega} - s_1 \Im\hat{\eta}_{\omega}) \int_{\partial D} \boldsymbol{\varphi}_1 =0, \medskip \\ \mathbf{S}_D[\boldsymbol{\varphi}_2] + (s_1 \Re\hat{\eta}_{\omega} + \Im\hat{\eta}_{\omega} ) \int_{\partial D} \boldsymbol{\varphi}_1 =0. \end{cases}\] Then we consider the following two systems with \(s_1\neq 0\), \[\begin{cases} \mathbf{S}_D[s_1 \boldsymbol{\varphi}_1] + (s_1\Re\hat{\eta}_{\omega} - s^2_1 \Im\hat{\eta}_{\omega}) \int_{\partial D} \boldsymbol{\varphi}_1 =0, \medskip \\ \int_{\partial D} s_1\boldsymbol{\varphi}_1(y) = s_1\mathbf{a}, \end{cases}\] and \[\begin{cases} \mathbf{S}_D[\boldsymbol{\varphi}_2] + (s_1 \Re\hat{\eta}_{\omega} + \Im\hat{\eta}_{\omega} ) \int_{\partial D} \boldsymbol{\varphi}_1 =0, \medskip \\ \int_{\partial D} \boldsymbol{\varphi}_2(y) = s_1\mathbf{a}. \end{cases}\] From Lemma 4, we have that the last two systems have the same solution. Thus, we obtain that \(s_1^2 = -1\). This is a contradiction with the fact that \(s_1\in \mathbb{R}\).

   Now we consider the case for \(s_1 = 0\). If \(s_1 = 0\), then we construct a new function \[\boldsymbol{\varphi}_3 = \boldsymbol{\varphi}_1 - \frac{\Re \hat{\eta}_{\omega}}{\Im \hat{\eta}_{\omega}} \boldsymbol{\varphi}_2.\] We can easily verify that \(\mathbf{S}_D[\boldsymbol{\varphi}_3]=0\) and \(\int_{\partial D} \boldsymbol{\varphi}_3 =\mathbf{a}\neq 0\). This is a contradiction with our assumption that the dimension of the kernel of the operator \(\mathbf{S}_D\) is \(1\). Finally, we have that \(\boldsymbol{\varphi}\) is a real function. Again, from the choice of \(\omega\), we obtain that \(\int_{\partial D} \boldsymbol{\varphi}=0\). Then we further conclude that \(\boldsymbol{\varphi}=0\).

\(\bullet\) Next, we consider the case that the operator \(\mathbf{S}_D\) is invertible. From 18 , there holds that \[\label{eq:scon} \mathbf{S}_D[\boldsymbol{\varphi}] =\mathbf{a}_1,\tag{22}\] where \(\mathbf{a}_1=(\mathbf{a}_{1,1},\mathbf{a}_{1,2})^t= -\hat{\eta}_{\omega} \int_{\partial D} \boldsymbol{\varphi}\). Let \(\mathbf{e}_1=(1,0)^t\) and \(\mathbf{e}_2=(0,1)^t\). Then, we define two functions \[\boldsymbol{\psi}_1=\mathbf{S}_D^{-1}[\mathbf{e}_1], \quad \boldsymbol{\psi}_2=\mathbf{S}_D^{-1}[\mathbf{e}_2],\] since the operator \(\mathbf{S}_D\) is invertible. It follows from 22 that the function \(\boldsymbol{\varphi}\) can be written as \[\boldsymbol{\varphi}=\mathbf{a}_{1,1}\boldsymbol{\psi}_1 + \mathbf{a}_{1,2}\boldsymbol{\psi}_2.\] Substituting the last equation into 18 gives that \[\label{eq:pin1} \mathbf{a}_{1,1}\big(\mathbf{e}_1 + \hat{\eta}_{\omega} \int_{\partial D}\boldsymbol{\psi}_1 \big) + \mathbf{a}_{1,2} \big(\mathbf{e}_2 + \hat{\eta}_{\omega} \int_{\partial D}\boldsymbol{\psi}_2\big)=0.\tag{23}\] We rewrite the identity 23 into the following matrix form \[\widetilde{B}\mathbf{a}_{1}=\mathbf{0},\] where \[\label{eq:detb} \widetilde{B} = \begin{pmatrix} 1 + \hat{\eta}_{\omega} \tilde{b}_{11} & \hat{\eta}_{\omega} \tilde{b}_{21}\\ \hat{\eta}_{\omega} \tilde{b}_{12}& 1 + \hat{\eta}_{\omega} \tilde{b}_{22} \end{pmatrix}, \quad \tilde{b}_{ij} = \int_{\partial D}\boldsymbol{\psi}_{i,j}, \;\; i,j=1,2.\tag{24}\] Direct calculation shows that \[\Im\det(\widetilde{B})= \Im{\hat{\eta}_{\omega}} \big( (\tilde{b}_{11} + \tilde{b}_{22} ) + 2 \Re{\hat{\eta}_{\omega}} (\tilde{b}_{11}\tilde{b}_{22} - \tilde{b}_{12}\tilde{b}_{21}) \big).\] Then, from the choice of \(\omega\), i.e., \(\Im{\hat{\eta}_{\omega}}\neq 0\) and \(\Re{\hat{\eta}_{\omega}} \neq -\frac{1}{2}\frac{\tilde{b}_{11} + \tilde{b}_{22}}{\tilde{b}_{11}\tilde{b}_{22} - \tilde{b}_{12}\tilde{b}_{21}}\), we have that \(\det(\widetilde{B})\neq 0\). Thus, we have \(\mathbf{a}_1=0\). Finally, we can directly conclude \(\boldsymbol{\varphi}=0\) from 22 due to the invertibility of the operator \(\mathbf{S}_D\).

The proof is completed. ◻

Proof. By the definition of the operator \(\mathbf{K}_{D,1}^{(1)}\) in Lemma 1, Green’s identity and the fact \(\partial_{\boldsymbol{\nu}_{\mathbf{x}}}\tilde{\boldsymbol{\vartheta}_i}=0\), we have that \[\begin{align} \int_{\partial D} \mathbf{K}_{D,1}^{(1)}[\boldsymbol{\varphi}](\mathbf{x})\cdot {\boldsymbol{\vartheta}_i}(\mathbf{x})\;\mathrm{d}s(\mathbf{x})&=\int_{\partial D}\int_{\partial D} \partial_{\boldsymbol{\nu}_{\mathbf{x}}}\mathbf{\Gamma}^{1,1}(\mathbf{x}-\mathbf{y})\boldsymbol{\varphi}(\mathbf{y})\;\mathrm{d}s(\mathbf{y})\cdot {\boldsymbol{\vartheta}_i}(\mathbf{x})\;\mathrm{d}s(\mathbf{x})\\ &= \int_{\partial D}\int_{ D} \mathcal{L}_{\lambda,\mu}\mathbf{\Gamma}^{1,1}(\mathbf{x}-\mathbf{y})\cdot \tilde{\boldsymbol{\vartheta}_i} (\mathbf{x})\;\mathrm{d}\mathbf{x}\cdot \boldsymbol{\varphi}(\mathbf{y})\;\mathrm{d}s(\mathbf{y}) \\ &=\tilde{\alpha}\rho \int_D\tilde{\boldsymbol{\vartheta}_i}(\mathbf{x}) \;\mathrm{d}\mathbf{x}\cdot \int_{\partial D}\boldsymbol{\varphi}(\mathbf{y})\;\mathrm{d}s(\mathbf{y}). \end{align}\] The proof is completed. ◻

Proof. By the definition of the operator \(\mathbf{K}_{D,1}^{(2)}\) in Lemma 1, Green’s identity and the fact \(\partial_{\boldsymbol{\nu}_{\mathbf{x}}}\tilde{\boldsymbol{\vartheta}_i}=0\), we have that \[\begin{align} \int_{\partial D} \mathbf{K}_{D,1}^{(2)}[\boldsymbol{\varphi}](\mathbf{x})\cdot {\boldsymbol{\vartheta}_i}(\mathbf{x})\;\mathrm{d}s(\mathbf{x})&=\int_{\partial D}\int_{\partial D} \partial_{\boldsymbol{\nu}_{\mathbf{x}}}\mathbf{\Gamma}^{1,2}(\mathbf{x}-\mathbf{y})\boldsymbol{\varphi}(\mathbf{y})\;\mathrm{d}s(\mathbf{y})\cdot {\boldsymbol{\vartheta}_i}(\mathbf{x})\;\mathrm{d}s(\mathbf{x})\\ &= \int_{\partial D}\int_{ D} \mathcal{L}_{\lambda,\mu}\mathbf{\Gamma}^{1,2}(\mathbf{x}-\mathbf{y})\cdot \tilde{\boldsymbol{\vartheta}_i} (\mathbf{x})\;\mathrm{d}\mathbf{x}\cdot \boldsymbol{\varphi}(\mathbf{y})\;\mathrm{d}s(\mathbf{y}) \\ &= -\rho\int_D{ \hat{\mathbf{S}}_D^\omega[\boldsymbol{\varphi}] \cdot \tilde{\boldsymbol{\vartheta}_i}} - \tilde{\alpha}\rho \ln\omega \int_D\tilde{\boldsymbol{\vartheta}_i} \cdot \int_{\partial D}\boldsymbol{\varphi}. \end{align}\] The proof is completed. ◻

Proof. With the help of the symmetry of the domain \(D\) and the change of variable, we have that \[\begin{align} \boldsymbol{\vartheta}_1 = \hat{\mathbf{S}}_D^\omega[\boldsymbol{\zeta}_1](\mathbf{x}) & = \int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\mathbf{x}-\mathbf{y})\boldsymbol{\zeta}_1(\mathbf{y})ds(\mathbf{y}) \\ & \overset{\mathbf{y}\to -\mathbf{y}}{=} \int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\mathbf{x}+\mathbf{y})\boldsymbol{\zeta}_1(-\mathbf{y})ds(\mathbf{y}) \\ & \overset{\mathbf{x}=-\widetilde{\mathbf{x}}}{=} \int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\widetilde{\mathbf{x}}-\mathbf{y})\boldsymbol{\zeta}_1(-\mathbf{y})ds(\mathbf{y}). \end{align}\] Thus, we can conclude that \(\boldsymbol{\zeta}_4(\mathbf{y}) = \boldsymbol{\zeta}_1(-\mathbf{y})\). Following the same argument, we can obtain that \(\boldsymbol{\zeta}_5(\mathbf{y}) = \boldsymbol{\zeta}_2(-\mathbf{y})\) and \(\boldsymbol{\zeta}_6(\mathbf{y}) = -\boldsymbol{\zeta}_3(-\mathbf{y})\). Then, direct calculation shows that \[\boldsymbol{\alpha}_{11}=\int_{\partial D_1} \boldsymbol{\zeta}_1(\mathbf{y})\;\mathrm{d}s(\mathbf{y})=\int_{\partial D_2} \boldsymbol{\zeta}_1(-\mathbf{y})\;\mathrm{d}s(\mathbf{y})=\int_{\partial D_2} \boldsymbol{\zeta}_4(\mathbf{y})\;\mathrm{d}s(\mathbf{y}) =\boldsymbol{\alpha}_{42}.\] Similarly, we can obtain the other equalities in ?? . By the definition of \(\boldsymbol{\alpha}_{i}\) in 27 , we get ?? . ◻

Proof. Since the domain \(D\) is symmetric with respect to the \(\mathbf{x}_1\)-axis, by the change of variables, we have that \[\label{eq:symx1} \begin{align} \boldsymbol{\vartheta}_2 = \hat{\mathbf{S}}_D^\omega[\boldsymbol{\zeta}_2](\mathbf{x}) & = \int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\mathbf{x}_1-\mathbf{y}_1,\mathbf{x}_2-\mathbf{y}_2)\boldsymbol{\zeta}_2(\mathbf{y}_1, \mathbf{y}_2)ds(\mathbf{y}) \\ & \overset{\mathbf{y}_2\to -\mathbf{y}_2}{=} \int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\mathbf{x}_1-\mathbf{y}_1,\mathbf{x}_2+\mathbf{y}_2)\boldsymbol{\zeta}_2(\mathbf{y}_1, -\mathbf{y}_2)ds(\mathbf{y}) \\ & \overset{\mathbf{x}_2=-\widetilde{\mathbf{x}}_2}{=} \int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\mathbf{x}_1-\mathbf{y}_1,-\widetilde{\mathbf{x}}_2+\mathbf{y}_2)\boldsymbol{\zeta}_2(\mathbf{y}_1, -\mathbf{y}_2)ds(\mathbf{y})\\ & = \int_{\partial D}(\widetilde{\mathbf{\Gamma}}^{0,1}+ \widetilde{\mathbf{\Gamma}}^{0,2})(\mathbf{x}_1-\mathbf{y}_1,\widetilde{\mathbf{x}}_2-\mathbf{y}_2)\widetilde{\boldsymbol{\zeta}}_2(\mathbf{y}_1, \mathbf{y}_2)ds(\mathbf{y}) \\ & = \int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\mathbf{x}_1-\mathbf{y}_1,\widetilde{\mathbf{x}}_2-\mathbf{y}_2)\widetilde{\boldsymbol{\zeta}}_2(\mathbf{y}_1, \mathbf{y}_2)ds(\mathbf{y}), \end{align}\tag{28}\] where \(\widetilde{\boldsymbol{\zeta}}_2\) is given in ?? and \[\left(\widetilde{\mathbf{\Gamma}}^{0,m}_{i, j}\right)_{i, j=1}^{2} = (-1)^{i}\left({\mathbf{\Gamma}}^{0,m}_{i, j}\right)_{i, j=1}^{2}, \quad m=1,2.\] The last identity in 28 follows from the fact that the first component of \(\boldsymbol{\vartheta}_2\) is zero. Thus, we have that \(\boldsymbol{\zeta}_5 = \widetilde{\boldsymbol{\zeta}}_2\). Following the same argument, we obtain that \[\boldsymbol{\zeta}_6 = \widetilde{\boldsymbol{\zeta}}_3, \quad \boldsymbol{\zeta}_4 = -\widetilde{\boldsymbol{\zeta}}_1,\] with \(\widetilde{\boldsymbol{\zeta}}_3\) and \(\widetilde{\boldsymbol{\zeta}}_1\) given in ?? . ◻

Proof. It is noted that the matrix \(A\) is symmetric, which follows from \[A_{ij}=\int_{\partial D} {\boldsymbol{\boldsymbol{\vartheta}}}_i \cdot \boldsymbol{\zeta}_j = \int_{\partial D} \hat{\mathbf{S}}_D^\omega [\boldsymbol{\zeta}_i] \cdot \boldsymbol{\zeta}_j =\int_{\partial D} \boldsymbol{\zeta}_i \cdot \hat{\mathbf{S}}_D^\omega [\boldsymbol{\zeta}_j] = \int_{\partial D} {\boldsymbol{\boldsymbol{\vartheta}}}_j \cdot \boldsymbol{\zeta}_i=A_{ji}.\] From Lemma 9, direct calculation shows that \[A_{ij}=A_{i+3,j+3}, \quad 1\leq i,j\leq 2, \quad A_{33}=A_{66},\] \[A_{ij}=A_{(i+3)(j-3)}, \quad i=1,2,\; j=4,5, \quad A_{36}=A_{63},\] \[A_{i3}= - A_{(i+3)6}, \quad A_{3i}= - A_{6(i+3)}, \quad i=1,2,\] \[A_{i6}= - A_{(i+3)3}, \quad A_{6i}= - A_{3(i+3)}, \quad i=1,2.\] By Lemma 10, we have that \[A_{ij}=A_{i+3,j+3}, \quad 2\leq i,j\leq 3, \quad A_{11}=A_{44},\] \[A_{ij}=A_{(i+3)(j-3)}, \quad i=2,3,\; j=5,6, \quad A_{14}=A_{41},\] \[A_{1i}= - A_{4(i+3)}, \quad A_{i1}= - A_{(i+3)4}, \quad i=2,3,\] \[A_{1i}= - A_{4(i-3)}, \quad A_{i1}= - A_{(i-3)4}, \quad i=5,6.\] From the above identities, we conclude ?? . The other identities directly follow from the values of \(\boldsymbol{\vartheta}_i\) and the expressions of the parameters defined in 27 . ◻

Proof. From the definition of \(\boldsymbol{\zeta}_i\) in 26 and the expression of the operator \(\hat{\mathbf{S}}_D^\omega\) in Lemma 1, we have that \[\label{eq:pral1} \hat{\mathbf{S}}_D^\omega[\boldsymbol{\zeta}_1] = (s_3 \ln\omega+s_4) \int_{\partial D} \boldsymbol{\zeta}_1+ \mathbf{S}_D[\boldsymbol{\zeta}_1] = \boldsymbol{\vartheta}_1.\tag{29}\] Next, we show that the function \(\boldsymbol{\zeta}_1\) has the following asymptotic expansion \[\label{eq:pral2} \boldsymbol{\zeta}_1(\mathbf{y}) = \boldsymbol{\zeta}_{1,0}(\mathbf{y}) + \sum_{j=1}^{\infty} \frac{1}{(\ln \omega)^j} \boldsymbol{\zeta}_{1,j}(\mathbf{y}).\tag{30}\] Substituting the last equation into 29 and comparing the order of the parameter \(\ln\omega\), we have that \[\label{eq:zeta1e1} \int_{\partial D}\boldsymbol{\zeta}_{1,0} = 0, \quad s_3\int_{\partial D}\boldsymbol{\zeta}_{1,1} + \mathbf{S}_D[\boldsymbol{\zeta}_{1,0}] = \boldsymbol{\vartheta}_1,\tag{31}\] and for \(j\geq 1\), \[\label{eq:zeta1e2} \begin{cases} s_3\int_{\partial D}\boldsymbol{\zeta}_{1,j+1} + s_4 \int_{\partial D}\boldsymbol{\zeta}_{1,j} + \mathbf{S}_D[\boldsymbol{\zeta}_{1,j}] = 0,\\ \int_{\partial D}\boldsymbol{\zeta}_{1,j} = \mathbf{t}_j, \end{cases}\tag{32}\] where \(\mathbf{t}_j\) are some vector-valued constants.

The existence of the solution \(\boldsymbol{\zeta}_{1,0}\) in 31 is guaranteed by Lemma 4. Next, we consider the system 32 . By Lemma 4, the solution \(\boldsymbol{\zeta}_{1,1}\) exists and satisfies the following estimate \[\|\boldsymbol{\zeta}_{1,1}\|_{L^2(\partial D)} + \left\vert s_3 \int_{\partial D}\boldsymbol{\zeta}_{1,2} + s_4 \mathbf{t}_1 \right\vert \leq C_{\mathcal{T}} \left\vert\mathbf{t}_1\right\vert,\] where \(C_{\mathcal{T}}\) is the norm of the operator \(\mathcal{T}\) defined in Lemma 4. From the last inequality and the triangle inequality, we have that \[\|\boldsymbol{\zeta}_{1,1}\|_{L^2(\partial D)} \leq C_{\mathcal{T}}\left\vert\mathbf{t}_1\right\vert \quad and \quad \left\vert\int_{\partial D}\boldsymbol{\zeta}_{1,2}\right\vert \leq \frac{C_{\mathcal{T}} + \left\vert s_4\right\vert}{\left\vert s_3\right\vert} \left\vert\mathbf{t}_1\right\vert.\] By the same argument, we can show that \[\|\boldsymbol{\zeta}_{1,j}\|_{L^2(\partial D)} \leq C_{\mathcal{T}} \left( \frac{C_{\mathcal{T}} + \left\vert s_4\right\vert}{\left\vert s_3\right\vert}\right)^{j-1}\left\vert\mathbf{t}_1\right\vert \quad and \quad \left\vert\int_{\partial D}\boldsymbol{\zeta}_{1,j+1}\right\vert \leq \left( \frac{C_{\mathcal{T}} + \left\vert s_4\right\vert}{\left\vert s_3\right\vert}\right)^{j} \left\vert\mathbf{t}_1\right\vert.\] Thus, if \(\frac{C_{\mathcal{T}} + \left\vert s_4\right\vert}{\left\vert s_3\right\vert|\ln\omega|} <1\), we conclude that the asymptotic expansion 30 converges with respect to the norm \(\|\cdot\|_{L^2(\partial D)}\). Therefore, we have \(\boldsymbol{\alpha}_1 \leq \mathcal{O}(1/\ln\omega)\) and the estimate of \(\boldsymbol{\alpha}_2\) and \(\boldsymbol{\alpha}_3\) follows from the same argument.

Next, we give the estimate in ?? . If the domain \(D\) is symmetric with respect to the origin, by Lemma 9 we have that similar to 31 , \[\label{eq:zeta1e3} \int_{\partial D}\boldsymbol{\zeta}_{4,0} = 0, \quad s_3\mathbf{t}_1 + \mathbf{S}_D[\boldsymbol{\zeta}_{4,0}] = \boldsymbol{\vartheta}_4,\tag{33}\] where \(\mathbf{t}_1=\int_{\partial D}\boldsymbol{\zeta}_{1,1}\). Adding 33 and 31 gives that \[\int_{\partial D}\boldsymbol{\zeta}_{1,0} + \boldsymbol{\zeta}_{4,0} = 0, \quad \mathbf{S}_D[\boldsymbol{\zeta}_{1,0} + \boldsymbol{\zeta}_{4,0}] + 2s_3 \mathbf{t}_1 - \begin{pmatrix} 1\\ 1 \end{pmatrix} = 0.\] By Lemma 4, we conclude \(\mathbf{t}_1=\frac{1}{2s_3} (1,1)^t\) since the operator \(\mathcal{T}\) is invertible. Finally, from the asymptotic expansion 30 , we obtain the estimate of \({\boldsymbol{\alpha}_1}\) in ?? . Following the same argument, we can obtain the estimate of \({\boldsymbol{\alpha}_2}\). The proof is completed. ◻

Proof. From the expression of the functions \({\boldsymbol{\vartheta}}_i\) and the definition of the entry \(B_{ij}\), we have ?? and \[B_{12}=B_{21}=B_{45}=B_{54}=0.\] Since the domain \(D\) is symmetric with respect to the origin and the \(\mathbf{x}_1\)-axis, we have that \[B_{23}=B_{32}=B_{56}=B_{65}=0.\] Finally, the direct calculation shows ?? . ◻

Proof. From the discussion above the theorem, we conclude that the system 4 shares resonant frequencies; that is, there exist resonant frequencies \(\omega\in\mathbb{C}\) such that \(\mathcal{A}(\omega,\delta)\) defined in 9 has a nontrivial kernel. Assume that \((\boldsymbol{\varphi}, \boldsymbol{\psi})^t\) is in the kernel, namely, \[\mathcal{A}(\omega,\delta)\begin{bmatrix} \boldsymbol{\varphi}\\ \boldsymbol{\psi} \end{bmatrix}=0.\] From the expression of the operator \(\mathcal{A}(\omega,\delta)\), the last equation is equivalent to \[\begin{cases} \mathbf{S}_D^{\tau \omega}[\boldsymbol{\varphi}]-\mathbf{S}_D^\omega [\boldsymbol{\psi}]=0,\\ \left(-\frac{1}{2}{\boldsymbol{I}}+\mathbf{K}_D^{\tau \omega,*}\right)[\boldsymbol{\varphi}]-\delta\left( \frac{1}{2}{\boldsymbol{I}}+\mathbf{K}^{\omega,*}_D\right)[\boldsymbol{\psi}]=0. \end{cases}\] By the asymptotic expansion of the operators \(\mathbf{S}_D^{\omega}\) and \(\mathbf{K}_D^{\tau \omega,*}\) in Lemma 1, we have that \[\label{eq:eee1} \hat{\mathbf{S}}_D^{\tau \omega} [\boldsymbol{\varphi}]-\hat{\mathbf{S}}^\omega [\boldsymbol{\psi}] =\mathcal{O}(\omega^2\ln{\omega}),\tag{37}\] and \[\label{eq:eee3} \begin{align} &\left(-\frac{1}{2}{\boldsymbol{I}} +\mathbf{K}_D^*+(\tau \omega)^2\ln{(\tau \omega)}\mathbf{K}_{D,1}^{(1)} + (\tau \omega)^2 \mathbf{K}_{D,1}^{(2)}\right)[\boldsymbol{\varphi}]-\delta \left( \frac{1}{2}{\boldsymbol{I}}+\mathbf{K}_D^*\right)[\boldsymbol{\psi}] \\ = & \mathcal{O}\left((\tau \omega)^4 \ln{(\tau \omega)}+\delta \omega^2\ln{\omega}\right). \end{align}\tag{38}\] From 38 , we have that \[\boldsymbol{\varphi}-\mathcal{O}((\tau \omega)^2\ln{(\tau\omega)}+\delta) \in \ker{\left(-\frac{1}{2}{\boldsymbol{I}} +\mathbf{K}_D^*\right)}.\] By Lemma 6, the function \(\boldsymbol{\varphi}\) can be written as \[\boldsymbol{\varphi}=\sum_{j=1}^6 d_j\boldsymbol{\xi}_j+\mathcal{O}((\tau \omega)^2\ln{(\tau\omega)}+\delta),\] where \(d_j\), \(1\leq j\leq 6\), are some constants needed to be determined. Then, from 26 and 37 , the function \(\boldsymbol{\psi}\) can be expressed by \[\boldsymbol{\psi}=\sum_{j=1}^6d_j\boldsymbol{\zeta}_j +\mathcal{O}( \omega^2\ln\omega+\delta).\] Substituting the last two equations into 38 gives that \[\label{eq:gere1} \begin{align} & \left(-\frac{1}{2}{\boldsymbol{I}} +\mathbf{K}_D^*\right)[\boldsymbol{\varphi}]+ (\tau\omega)^2\ln{(\tau\omega)}\mathbf{K}_{D,1}^{(1)}\left(\sum_{j=1}^6 d_j\boldsymbol{\xi}_j \right) + (\tau\omega)^2 \mathbf{K}_{D,1}^{(2)}\left(\sum_{j=1}^6 d_j\boldsymbol{\xi}_j \right) \\ &-\delta \left(\frac{1}{2}{\boldsymbol{I}}+\mathbf{K}_D^*\right)\left( \sum_{j=1}^6 d_j\boldsymbol{\zeta}_j \right) =\mathcal{O}\left( (\tau \omega)^4 (\ln{(\tau \omega)})^2 +\delta \omega^2 \ln{\omega} +\delta^2 \right). \end{align}\tag{39}\] Multiplying \({\boldsymbol{\vartheta}}_i\), \(1\leq i\leq 6\), on both sides of 39 , integrating on \(\partial D\) and with the help of Lemma 5 gives that \[\begin{align} &(\tau \omega)^2\ln{(\tau\omega)} \int_{\partial D} \mathbf{K}_{D,1}^{(1)}\left(\sum_{j=1}^6 d_j\boldsymbol{\xi}_j\right) \cdot {\boldsymbol{\boldsymbol{\vartheta}}}_i\; + (\tau \omega)^2 \int_{\partial D} \mathbf{K}_{D,1}^{(2)}\left(\sum_{j=1}^6 d_j\boldsymbol{\xi}_j\right) \cdot {\boldsymbol{\boldsymbol{\vartheta}}}_i\; \\ & -\delta \int_{\partial D} \left(\frac{1}{2}{\boldsymbol{I}}+\mathbf{K}_D^*\right)\left( \sum_{j=1}^6 d_j\boldsymbol{\zeta}_j \right) \cdot{\boldsymbol{\vartheta}}_i\; =\mathcal{O}\left( (\tau \omega)^4 (\ln{(\tau \omega)})^2 +\delta \omega^2 \ln{\omega} +\delta^2 \right). \end{align}\] By Lemmas 5, 7 and 8, the last equation can be written as, for \(1\leq i\leq 6\), \[\begin{align} & \sum_{j=1}^6 d_j\left( \rho (\tau \omega)^2 \int_D\tilde{\boldsymbol{\vartheta}_i}(\mathbf{y}) \cdot \tilde{\boldsymbol{\vartheta}_j}(\mathbf{y}) + \delta \int_{\partial D} {\boldsymbol{\boldsymbol{\vartheta}}}_i(\mathbf{y}) \cdot \boldsymbol{\zeta}_j(\mathbf{y}) \right)\\ &= \mathcal{O}\left( (\tau \omega)^4 (\ln{(\tau \omega)})^2 +\delta \omega^2 \ln{\omega} +\delta^2 \right). \end{align}\] Using the notations given in Proposition 1 and Lemma 12, namely, \(A_{ij}=\int_{\partial D} {\boldsymbol{\boldsymbol{\vartheta}}}_i(\mathbf{y}) \cdot \boldsymbol{\zeta}_j(\mathbf{y}) \;\mathrm{d}s(\mathbf{y})\) and \(B_{ij}=\int_D\tilde{\boldsymbol{\vartheta}_i}(\mathbf{y}) \cdot \tilde{\boldsymbol{\vartheta}_j}(\mathbf{y}) \;\mathrm{d}\mathbf{y}\), the leading terms of the last equations can be written as \[\label{eq:adb1} \left(\rho (\tau \omega)^2 B + \delta A \right)\mathbf{d}= 0,\tag{40}\] where \(\mathbf{d}= (d_1, d_2, d_3, d_4, d_5, d_6)^t\). Then, by Lemma 9, Proposition 1, and Lemma 12, 40 can be decomposed into two sub-equations \[\label{eq:adb2} \left(\rho (\tau \omega)^2 B_i + \delta A_i \right)\tilde{\mathbf{d}}_i = 0, \quad i=1,2,\tag{41}\] where \(A_1\) and \(B_1\) are 2-by-2 matrices, and \(A_2\) and \(B_2\) are 4-by-4 matrices.

For matrices \(A_1\) and \(B_1\), they are given by \[B_1=\begin{pmatrix} |D_1|& 0\\ 0& |D_1| \end{pmatrix}, \quad \quad A_1=\begin{pmatrix} \boldsymbol{\alpha}_{21,2} & \boldsymbol{\alpha}_{22,2}\\ \boldsymbol{\alpha}_{22,2} & \boldsymbol{\alpha}_{21,2} \end{pmatrix}.\] Direct calculation shows that the eigensystem of the matrix \(A_1\) is given by \[\boldsymbol{\alpha}_{21,2}-\boldsymbol{\alpha}_{22,2}, \; (1, -1)^t; \quad \boldsymbol{\alpha}_{21,2}+\boldsymbol{\alpha}_{22,2}, \; (1, 1)^t.\] Thus, from equation 41 , the first resonant frequency satisfies \[\omega_1 = \frac{1}{\tau} \sqrt{\frac{\delta (\boldsymbol{\alpha}_{22,2}-\boldsymbol{\alpha}_{21,2})}{\rho |D_1|}}(1 + o(1)),\] The corresponding eigenfunctions are given by \[\boldsymbol{\varphi}_1 =\boldsymbol{\xi}_2 -\boldsymbol{\xi}_5+\mathcal{O}((\tau \omega_1)^2\ln{(\tau\omega_1)}), \quad \boldsymbol{\psi}_1 =\boldsymbol{\zeta}_2 -\boldsymbol{\zeta}_5+\mathcal{O}( \omega_1^2\ln\omega_1),\] where we used the fact \(\delta = \mathcal{O}((\tau \omega_1)^2)\). From Lemma 11, the leading order of the second resonant frequency satisfies \[2 s_3 \rho |D_1| (\tau \omega_{2,1})^2 \ln{ \omega_{2,1}} + \delta = 0.\] The corresponding eigenfunctions are given by \[\boldsymbol{\varphi}_2 =\boldsymbol{\xi}_2 +\boldsymbol{\xi}_5+\mathcal{O}((\tau \omega_2)^2\ln{(\tau\omega_2)}), \quad \boldsymbol{\psi}_2 =\boldsymbol{\zeta}_2 +\boldsymbol{\zeta}_5+\mathcal{O}( \omega_2^2\ln\omega_2),\] which follows from the fact \(\delta = \mathcal{O}((\tau \omega_2)^2\ln \omega_2)\).

For matrices \(A_2\) and \(B_2\), from Lemma 9, they have the following expressions: \[B_2 =\begin{pmatrix} |D_1| & \gamma_1 & 0& 0\\ \gamma_1 & \gamma_2 & 0 & 0\\ 0&0&|D_1|&-\gamma_1&\\ 0&0& -\gamma_1 &\gamma_2 \end{pmatrix}, \quad A_2=\begin{pmatrix} \boldsymbol{\alpha}_{11,1} &\boldsymbol{\alpha}_{31,1} &\boldsymbol{\alpha}_{12,1} &-\boldsymbol{\alpha}_{32,1} \\ \boldsymbol{\alpha}_{31,1} &\beta_{31} &\boldsymbol{\alpha}_{32,1} &\beta_{32}\\ \boldsymbol{\alpha}_{12,1} &\boldsymbol{\alpha}_{32,1} &\boldsymbol{\alpha}_{11,1} &-\boldsymbol{\alpha}_{31,1} \\ -\boldsymbol{\alpha}_{32,1} &\beta_{32} &-\boldsymbol{\alpha}_{31,1} &\beta_{31} \end{pmatrix}.\] To find the resonant frequencies, we consider the following equation instead of 41 : \[\label{eq:adb3} \left(\rho (\tau \omega)^2 + \delta B_2^{-1} A_2 \right)\tilde{\mathbf{d}}_2 = 0.\tag{42}\] The first eigenvalue of the matrix \(B_2^{-1} A_2\) is given by \[\iota_1 = \frac{\left((\beta_{32} - \beta_{31}) |D_1| - (\boldsymbol{\alpha}_{11,1} + \boldsymbol{\alpha}_{12,1}) \gamma_{2} + 2 (\boldsymbol{\alpha}_{31,1} + \boldsymbol{\alpha}_{32,1} ) \gamma_{1} + \sqrt{\iota_{1,1}} \right)}{-2\left( |D_1| \gamma_{2} - \gamma_{1}^{2} \right) } ,\] where \[\label{eq:iota11} \begin{align} \iota_{1,1} = &((\beta_{31} - \beta_{32}) |D_1| + (\boldsymbol{\alpha}_{11,1} + \boldsymbol{\alpha}_{12,1}) \gamma_{2} - 2(\boldsymbol{\alpha}_{31,1} + \boldsymbol{\alpha}_{32,1} ) \gamma_{1})^2 + \\ & 4 ((\boldsymbol{\alpha}_{31,1} + \boldsymbol{\alpha}_{32,1})^2 - (\boldsymbol{\alpha}_{11,1} + \boldsymbol{\alpha}_{12,1})(\beta_{31} - \beta_{32})) (|D_1|\gamma_{2} - \gamma_{1}^2). \end{align}\tag{43}\] By Lemma 11, the eigenvalue \(\iota_1\) has the following asymptotic expansion \[\iota_1 = \frac{ (\beta_{32} - \beta_{31}) |D_1|}{-\left( |D_1| \gamma_{2} - \gamma_{1}^{2} \right) } \big(1 + o(1) \big).\] Thus, from 42 , the third resonant frequency satisfies \[\omega_3 = \frac{1}{\tau} \sqrt{\frac{ (\beta_{32} - \beta_{31}) |D_1|}{\left( |D_1| \gamma_{2} - \gamma_{1}^{2} \right)\rho } \delta }(1 + o(1)).\] The corresponding eigenfunctions are given by \[\begin{align} \boldsymbol{\varphi}_3 &= e_{3,1}\boldsymbol{\xi}_1 + e_{3,3}\boldsymbol{\xi}_3 + e_{3,4}\boldsymbol{\xi}_4 + e_{3,6}\boldsymbol{\xi}_6 + \mathcal{O}((\tau \omega_3)^2\ln{(\tau\omega_3)}), \\ \boldsymbol{\psi}_3 &= e_{3,1}\boldsymbol{\zeta}_1 + e_{3,3}\boldsymbol{\zeta}_3 + e_{3,4}\boldsymbol{\zeta}_4 + e_{3,6}\boldsymbol{\zeta}_6 + \mathcal{O}( \omega_3^2\ln\omega_3), \end{align}\] where \[\label{def-e3} e_{3,1} = e_{3,4} =\frac{\gamma_1}{|D_1|}(1 + o(1)), \quad e_{3,3}=-1 +o(1), \quad e_{3,6}=1 +o(1).\tag{44}\]

The second eigenvalue of the matrix \(B_2^{-1} A_2\) is given by \[\iota_2 = \frac{\left((\beta_{32} - \beta_{31}) |D_1| - (\boldsymbol{\alpha}_{11,1} + \boldsymbol{\alpha}_{12,1}) \gamma_{2} + 2 (\boldsymbol{\alpha}_{31,1} + \boldsymbol{\alpha}_{32,1} ) \gamma_{1} - \sqrt{\iota_{1,1}} \right)}{-2\left( |D_1| \gamma_{2} - \gamma_{1}^{2} \right) } ,\] where \(\iota_{1,1}\) is defined in 43 . By Lemma 11 and direct calculation, the eigenvalue \(\iota_2\) has the following asymptotic expansion \[\iota_2 = \frac{\boldsymbol{\alpha}_{11,1} + \boldsymbol{\alpha}_{12,1}}{|D_1|} \big(1 + o(1) \big).\] Thus, from 42 and Lemma 11, the leading order of the fourth resonant frequency satisfies \[2 s_3 \rho |D_1| (\tau \omega_{4,1})^2 \ln{ \omega_{4,1}} + \delta = 0.\] The corresponding eigenfunctions are given by \[\begin{align} \boldsymbol{\varphi}_4 &= e_{4,1}\boldsymbol{\xi}_1 + e_{4,3}\boldsymbol{\xi}_3 + e_{4,4}\boldsymbol{\xi}_4 + e_{4,6}\boldsymbol{\xi}_6 + \mathcal{O}((\tau \omega_4)^2\ln{(\tau\omega_4)}), \\ \boldsymbol{\psi}_4 &= e_{4,1}\boldsymbol{\zeta}_1 + e_{4,3}\boldsymbol{\zeta}_3 + e_{4,4}\boldsymbol{\zeta}_4 + e_{4,6}\boldsymbol{\zeta}_6 + \mathcal{O}( \omega_4^2\ln\omega_4), \end{align}\] where \[\label{def-e4} e_{4,1} = e_{4,4} =\frac{\beta_{32} - \beta_{31}}{\boldsymbol{\alpha}_{31,1} - \boldsymbol{\alpha}_{32,1}}(1 + o(1)), \quad e_{4,3}=-1 +o(1), \quad e_{4,6}=1 +o(1).\tag{45}\]

The third eigenvalue of the matrix \(B_2^{-1} A_2\) is given by \[\iota_3 = \frac{(-\beta_{32} - \beta_{31}) |D_1| + (\boldsymbol{\alpha}_{12,1} - \boldsymbol{\alpha}_{11,1}) \gamma_{2} + 2(\boldsymbol{\alpha}_{31,1} - \boldsymbol{\alpha}_{32,1}) \gamma_{1} + \sqrt{\iota_{3,1}} }{-2\left( |D_1| \gamma_{2} - \gamma_{1}^{2} \right) } ,\] where \[\label{eq:iota31} \begin{align} \iota_{3,1} = &((\beta_{31} + \beta_{32}) |D_1| + (\boldsymbol{\alpha}_{11,1} - \boldsymbol{\alpha}_{12,1}) \gamma_{2} - 2(\boldsymbol{\alpha}_{31,1} - \boldsymbol{\alpha}_{32,1} ) \gamma_{1})^2 + \\ & 4 ((\boldsymbol{\alpha}_{31,1} - \boldsymbol{\alpha}_{32,1})^2 - (\boldsymbol{\alpha}_{11,1} - \boldsymbol{\alpha}_{12,1})(\beta_{31} + \beta_{32})) (|D_1|\gamma_{2} - \gamma_{1}^2). \end{align}\tag{46}\] By Proposition 2, the eigenvalue \(\iota_3\) has the following asymptotic expansion \[\iota_3 = \frac{ (\boldsymbol{\alpha}_{12,1} - \boldsymbol{\alpha}_{11,1}) \gamma_2 }{-(|D_{1}|\, \gamma_{2} - \gamma_{1}^2) } (1 + o(1)),\] Thus, from 42 , the fifth resonant frequency satisfies \[\omega_5 = \frac{1}{\tau} \sqrt{\frac{ (\boldsymbol{\alpha}_{12,1} - \boldsymbol{\alpha}_{11,1}) \gamma_2 }{(|D_{1}|\, \gamma_{2} - \gamma_{1}^2)\rho } \delta }(1 + o(1)).\] The corresponding eigenfunctions are given by \[\begin{align} \boldsymbol{\varphi}_5 &= e_{5,1}\boldsymbol{\xi}_1 + e_{5,3}\boldsymbol{\xi}_3 + e_{5,4}\boldsymbol{\xi}_4 + e_{5,6}\boldsymbol{\xi}_6 + \mathcal{O}((\tau \omega_5)^2\ln{(\tau\omega_5)}), \\ \boldsymbol{\psi}_5 &= e_{5,1}\boldsymbol{\zeta}_1 + e_{5,3}\boldsymbol{\zeta}_3 + e_{5,4}\boldsymbol{\zeta}_4 + e_{5,6}\boldsymbol{\zeta}_6 + \mathcal{O}( \omega_5^2\ln\omega_5), \end{align}\] where \[\label{def-e5} e_{5,1} = -e_{5,4} = 1 +o(1), \quad e_{5,3}= e_{5,6}= \frac{(\boldsymbol{\alpha}_{11,1} - \boldsymbol{\alpha}_{12,1})\gamma_{2}}{ -(\boldsymbol{\alpha}_{31,1} - \boldsymbol{\alpha}_{32,1})\gamma_{2} + (\beta_{31} + \beta_{32}) \gamma_{1} }(1 + o(1)).\tag{47}\]

The fourth eigenvalue of the matrix \(B_2^{-1} A_2\) is given by \[\iota_4 = \frac{(-\beta_{32} - \beta_{31}) |D_1| + (\boldsymbol{\alpha}_{12,1} - \boldsymbol{\alpha}_{11,1}) \gamma_{2} + 2(\boldsymbol{\alpha}_{31,1} - \boldsymbol{\alpha}_{32,1}) \gamma_{1} - \sqrt{\iota_{3,1}} }{-2\left( |D_1| \gamma_{2} - \gamma_{1}^{2} \right) } ,\] where \(\iota_{3,1}\) is defined in 46 . By Proposition 2, the eigenvalue \(\iota_4\) has the following asymptotic expansion \[\iota_4 = \frac{\beta_{31} + \beta_{32}}{\gamma_{2}} (1 + o(1)).\] Thus, from 42 , the sixth resonant frequency satisfies \[\omega_6 = \frac{1}{\tau} \sqrt{\frac{\beta_{31} + \beta_{32}}{-\gamma_{2} \rho} \delta }(1 + o(1)).\] The corresponding eigenfunctions are given by \[\begin{align} \boldsymbol{\varphi}_6 &= e_{6,1}\boldsymbol{\xi}_1 + e_{6,3}\boldsymbol{\xi}_3 + e_{6,4}\boldsymbol{\xi}_4 + e_{6,6}\boldsymbol{\xi}_6 + \mathcal{O}((\tau \omega_6)^2\ln{(\tau\omega_6)}), \\ \boldsymbol{\psi}_6 &= e_{6,1}\boldsymbol{\zeta}_1 + e_{6,3}\boldsymbol{\zeta}_3 + e_{6,4}\boldsymbol{\zeta}_4 + e_{6,6}\boldsymbol{\zeta}_6 + \mathcal{O}( \omega_6^2\ln\omega_6), \end{align}\] where \[\label{def-e6} e_{6,1} = -e_{6,4} =-\frac{ \gamma_2 }{ \gamma_1 }(1 + o(1)), \quad e_{6,3}= e_{6,6}=1 +o(1).\tag{48}\] The proof is completed. ◻

Proof. By using the Sobolev embedding theorem and classical \(W^{2,q}\)-estimate for elliptic systems, for \(q>2\), we obtain \[\begin{align} \|\nabla \mathbf{w}_{i,2}\|_{L^\infty(\Omega\setminus\Omega_{2R_0})}&\leq C\|\mathbf{w}_{i,2}\|_{W^{2,q}(\Omega\setminus\Omega_{2R_0})}\nonumber\\ &\leq C\big(\|\mathbf{w}_{i,2}\|_{L^{2}(\Omega\setminus\Omega_{R_0})}+\|{\boldsymbol{f}}_i\|_{L^{\infty}(\partial B_R)}\big). \end{align}\] This implies that we only need to consider \[\label{eq:bddu12} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}\mathbf{w}_{i,2}=0\quad &\text{in} \quad \Omega_{2R_0},\\ \mathbf{w}_{i,2}=0\quad &\text{on} \quad \Gamma_+^{2R_0}\cup\Gamma_-^{2R_0}, \end{cases}\tag{58}\] and estimate \(\nabla\mathbf{w}_{i,2}\) in the narrow region \(\Omega_{2R_0}\). Denote \[\Gamma_+^{2R_0}:=\left\{({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)\in\mathbb{R}^2: ~{\boldsymbol{x}}_2=\frac{\varepsilon}{2}+\mathcal{H}_{1}({\boldsymbol{x}}_1),~|{\boldsymbol{x}}_1|\leq 2R_0\right\}\] and \[\Gamma_-^{2R_0}:=\left\{({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)\in\mathbb{R}^2: ~{\boldsymbol{x}}_2=-\frac{\varepsilon}{2}+\mathcal{H}_{2}({\boldsymbol{x}}_1),~|{\boldsymbol{x}}_1|\leq 2R_0\right\}.\] Multiplying the equation in 58 by \(\mathbf{w}_{i,2}\) and integrating by parts, we have \[\begin{align} &\int_{\Omega_{2R_0}}\left(\mathbb{C}e(\mathbf{w}_{i,2}),e(\mathbf{w}_{i,2})\right)\\ &=\int_{\substack{|{\boldsymbol{x}}_1|=2R_0\\ -\frac{\varepsilon}{2}+\mathcal{H}_{2}({\boldsymbol{x}}_1)<{\bf x}_2<\frac{\varepsilon}{2}+\mathcal{H}_{1}({\boldsymbol{x}}_1)}}\mathbf{w}_{i,2}\big(\lambda(\nabla\cdot \mathbf{w}_{i,2})+\mu(\nabla \mathbf{w}_{i,2}+(\nabla\mathbf{w}_{i,2})^{t})\big)\frac{\mathbf{x}}{r}\\ &\leq C\int_{\substack{|{\boldsymbol{x}}_1|=2R_0\\ -\frac{\varepsilon}{2}+\mathcal{H}_{2}({\boldsymbol{x}}_1)<{\bf x}_2<\frac{\varepsilon}{2}+\mathcal{H}_{1}({\boldsymbol{x}}_1)}}(|\mathbf{w}_{i,2}|^2+|\nabla \mathbf{w}_{i,2}|^2). \end{align}\] By using [51], we have \[\begin{align} \int_{\Omega_{2R_0}}|\nabla\mathbf{w}_{i,2}|^2\leq C\int_{\Omega_{2R_0}}|e(\mathbf{w}_{i,2})|^2\leq C\int_{\Omega_{2R_0}}\left(\mathbb{C}e(\mathbf{w}_{i,2}),e(\mathbf{w}_{i,2})\right). \end{align}\] Thus, we obtain \[\begin{align} \int_{\Omega_{2R_0}}|\nabla\mathbf{w}_{i,2}|^2\leq C\int_{\substack{|{\boldsymbol{x}}_1|=2R_0\\ -\frac{\varepsilon}{2}+\mathcal{H}_{2}({\boldsymbol{x}}_1)<{\bf x}_2<\frac{\varepsilon}{2}+\mathcal{H}_{1}({\boldsymbol{x}}_1)}}(|\mathbf{w}_{i,2}|^2+|\nabla \mathbf{w}_{i,2}|^2). \end{align}\] Recalling that \(\mathbf{w}_{i,2}=0\) on \(\partial D_1\), we have, for any \(\mathbf{x}\in\Omega_{3R_0}\setminus\Omega_{3R_0/2}\), \[\begin{align} |\mathbf{w}_{i,2}({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)|&=|\mathbf{w}_{i,2}({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)-\mathbf{w}_{i,2}({\boldsymbol{x}}_1,\frac{\varepsilon}{2}+\mathcal{H}_{1}({\boldsymbol{x}}_1))|\\ &\leq C(\varepsilon+\kappa|{\boldsymbol{x}}_1|^2)\|\nabla \mathbf{w}_{i,2}\|_{L^\infty(\Omega_{3R_0}\setminus\Omega_{3R_0/2})}\\ &\leq C\big(\|\mathbf{w}_{i,2}\|_{L^{2}(\Omega\setminus\Omega_{R_0})}+\|{\boldsymbol{f}}_i\|_{L^{\infty}(\partial B_R)}\big), \end{align}\] where we used the Sobolev embedding theorem and classical \(W^{2,q}\)-estimate for elliptic systems in the last inequality. Hence, we derive \[\begin{align} \int_{\Omega_{2R_0}}|\nabla\mathbf{w}_{i,2}|^2&\leq C\int_{\substack{|{\boldsymbol{x}}_1|=2R_0\nonumber\\ -\frac{\varepsilon}{2}+\mathcal{H}_{2}({\boldsymbol{x}}_1)<{\bf x}_2<\frac{\varepsilon}{2}+\mathcal{H}_{1}({\boldsymbol{x}}_1)}}(|\mathbf{w}_{i,2}|^2+|\nabla \mathbf{w}_{i,2}|^2)\\ &\leq C\big(\|\mathbf{w}_{i,2}\|_{L^{2}(\Omega\setminus\Omega_{R_0})}^2+\|{\boldsymbol{f}}_i\|_{L^{\infty}(\partial B_R)}^2\big). \end{align}\] With this inequality in hand, by adapting the iterative arguments used in the proof of [35], we derive \[|\nabla\mathbf{w}_{i,2}(\mathbf{x})|\leq C\big(\|\mathbf{w}_{i,2}\|_{L^2(\Omega)}+\|{\boldsymbol{f}}_{i}\|_{L^{\infty}(\partial B_R)}\big),\quad\mathbf{x}\in\Omega_{2R_0}.\] The proof is finished. ◻

Proof of Proposition 2. Step 1: Proof of ?? . Define \[\mathcal{C}_{R_0}:=\left\{({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)\in \mathbb{R}^{2}: -\frac{\varepsilon}{2}+2\min_{|{\boldsymbol{x}}_1|=R_0}\mathcal{H}_{2}({\boldsymbol{x}}_1)<{\bf x}_2<\frac{\varepsilon}{2}+2\max_{|{\boldsymbol{x}}_1|=R_0}\mathcal{H}_{1}({\boldsymbol{x}}_1),~|{\boldsymbol{x}}_1|<R_0\right\}.\] Then combining with the gradient estimates for elliptic systems, we have \[\begin{align} \label{est-alp111} \boldsymbol{\alpha}_{11,1}=\int_{\partial D}\frac{\partial\mathbf{w}_1}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_1&=\int_{\partial D_1\cap\mathcal{C}_{R_0}}\frac{\partial\mathbf{w}_1}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_1+\int_{\partial D_1\setminus\mathcal{C}_{R_0}}\frac{\partial\mathbf{w}_1}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_1\nonumber\\ &=\int_{\partial D_1\cap\mathcal{C}_{R_0}}\frac{\partial\mathbf{w}_1}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_1+\mathcal{O}(1). \end{align}\tag{59}\] By using \[\begin{align} \label{def-nu-wi} \frac{\partial \mathbf{w}_{i}}{\partial\boldsymbol{\nu}}=\lambda(\nabla\cdot \mathbf{w}_{i}){\boldsymbol{n}}+\mu(\nabla \mathbf{w}_{i}+(\nabla\mathbf{w}_{i})^{t}){\boldsymbol{n}}, \end{align}\tag{60}\] where \[{\boldsymbol{n}}=(n_1,n_2),\quad n_1=\frac{\partial_{{\boldsymbol{x}}_1}\mathcal{H}_{1}({\boldsymbol{x}}_1)}{\sqrt{1+|\partial_{{\boldsymbol{x}}_1}\mathcal{H}_{1}({\boldsymbol{x}}_1)|^2}},\quad n_2=-\frac{1}{\sqrt{1+|\partial_{{\boldsymbol{x}}_1}\mathcal{H}_{1}({\boldsymbol{x}}_1)|^2}},\] we have \[\begin{align} \label{int-w1} &\int_{\partial D_1\cap\mathcal{C}_{R_0}}\frac{\partial\mathbf{w}_1}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_1\nonumber\\ &=\int_{\partial D_1\cap\mathcal{C}_{R_0}}\left(\lambda(\partial_{{\boldsymbol{x}}_1}w_{1}^1+\partial_{{\boldsymbol{x}}_2}w_{1}^2)n_1+\mu(2n_1\partial_{{\boldsymbol{x}}_1}w_{1}^1+n_2(\partial_{{\boldsymbol{x}}_1}w_{1}^2+\partial_{{\boldsymbol{x}}_2}w_{1}^1))\right). \end{align}\tag{61}\] In view of Proposition 3 and 57 , we find that the biggest term on the right-hand side of 61 is \[\mu\int_{\partial D_1\cap\mathcal{C}_{R_0}}n_2\partial_{{\boldsymbol{x}}_2}w_{1}^1,\] and the remaining terms are bounded by \(\mathcal{O}(1)\). By using a direct calculation together with 51 , we obtain \[\begin{align} \mu\int_{\partial D_1\cap\mathcal{C}_{R_0}}n_2\partial_{{\boldsymbol{x}}_2}w_{1}^1&=-\mu\int_{|{\boldsymbol{x}}_1|\leq R_0}\frac{1}{\delta({\boldsymbol{x}}_1)}\;d{\boldsymbol{x}}_1+\mathcal{O}(1)\\ &=-\mu\int_{|x_1|\leq R_0}\frac{1}{\varepsilon+\kappa|x_1|^2}\;dx_1+\varepsilon^{\frac{\alpha-1}{2}}\mathcal{O}(1)\\ &=-\frac{2\mu}{\sqrt{\kappa\varepsilon}}\int_{0}^{R_0\sqrt{\frac{\kappa}{\varepsilon}}}\frac{1}{1+r^2}\;dr+\varepsilon^{\frac{\alpha-1}{2}}\mathcal{O}(1)=-\frac{\mu\pi}{\sqrt{\kappa\varepsilon}}+\varepsilon^{\frac{\alpha-1}{2}}\mathcal{O}(1). \end{align}\] Thus, we obtain \[\boldsymbol{\alpha}_{11,1}=-\frac{\mu\pi}{\sqrt{\kappa\varepsilon}}+\varepsilon^{\frac{\alpha-1}{2}}\mathcal{O}(1).\] Similar to 59 , we derive \[\boldsymbol{\alpha}_{12,1}=\int_{\partial D}\frac{\partial\mathbf{w}_4}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_1=\int_{\partial D_1\cap\mathcal{C}_{R_0}}\frac{\partial\mathbf{w}_4}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_1+\mathcal{O}(1).\] By replicating the proof as above, we get that the leading order term of \(\boldsymbol{\alpha}_{12,1}\) is \[\begin{align} \mu\int_{\partial D_1\cap\mathcal{C}_{R_0}}n_2\partial_{{\boldsymbol{x}}_2}w_{4}^1&=\mu\int_{|{\boldsymbol{x}}_1|\leq R_0}\frac{1}{\delta({\boldsymbol{x}}_1)}\;d{\boldsymbol{x}}_1+\mathcal{O}(1)\\ &=\frac{\mu\pi}{\sqrt{\kappa\varepsilon}}+\varepsilon^{\frac{\alpha-1}{2}}\mathcal{O}(1). \end{align}\] This gives \[\boldsymbol{\alpha}_{12,1}=\frac{\mu\pi}{\sqrt{\kappa\varepsilon}}+\varepsilon^{\frac{\alpha-1}{2}}\mathcal{O}(1).\] The estimates of \(\boldsymbol{\alpha}_{21,2}\) and \(\boldsymbol{\alpha}_{22,2}\) are proved in a similar manner, we thus omit the details.

Step 2: Proof of ?? and ?? . Note that \[\begin{align} \beta_{3k}=\int_{\partial D_k}\frac{\partial\mathbf{w}_3}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_3=\int_{\partial D_k\cap\mathcal{C}_{R_0}}\frac{\partial\mathbf{w}_3}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\vartheta}}_3+\mathcal{O}(1),\quad k=1,2. \end{align}\] It follows from 57 that \[\begin{align} |\partial_{\mathbf{x}_1}\mathbf{v}_{3}^1|\leq C|\mathbf{x}_1|,\quad|\partial_{\mathbf{x}_2}\mathbf{v}_{3}^1|\leq C,\quad |\partial_{\mathbf{x}_1}\mathbf{v}_{3}^2|\leq C,\quad\partial_{x_2}\mathbf{v}_{3}^2=-\frac{\mathbf{x}_1}{\delta(\mathbf{x}_1)}. \end{align}\] Then by using Proposition 3 and 60 , we obtain \[\begin{align} \beta_{3k}&=\int_{\partial D_k\cap\mathcal{C}_{R_0}}\Big(\lambda(\partial_{{\boldsymbol{x}}_1}w_{3}^1+\partial_{{\boldsymbol{x}}_2}w_{3}^2)(n_1\mathbf{x}_2-n_2\mathbf{x}_1)\\ &\quad+\mu\big(\mathbf{x}_2(2n_1\partial_{{\boldsymbol{x}}_1}w_{3}^1+n_2(\partial_{{\boldsymbol{x}}_1}w_{3}^2+\partial_{{\boldsymbol{x}}_2}w_{1}^1))-\mathbf{x}_1(n_1(\partial_{{\boldsymbol{x}}_1}w_{3}^2+\partial_{{\boldsymbol{x}}_2}w_{1}^1)+2n_2\partial_{{\boldsymbol{x}}_2}w_{3}^2)\big)\Big)\\ &=\mathcal{O}(1),\quad k=1,2. \end{align}\] Similarly, we obtain \(\boldsymbol{\alpha}_{31,1}=\beta_{11}=\mathcal{O}(1)\) and \(\boldsymbol{\alpha}_{32,1}=-\beta_{12}=\mathcal{O}(1)\).

It follows from ?? and 27 that \[\begin{align} \label{beta3132} \beta_{31}+\beta_{32}=A_{33}+A_{36}=\int_{\partial D_1}\frac{\partial(\mathbf{w}_3+\mathbf{w}_6)}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}\boldsymbol{\varkappa}_3=\int_{\partial D_2}\frac{\partial(\mathbf{w}_3+\mathbf{w}_6)}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}\boldsymbol{\varkappa}_3. \end{align}\tag{62}\] By using 52 and \(\boldsymbol{\zeta}_6(-\mathbf{x})=-\boldsymbol{\zeta}_3(\mathbf{x})\), we obtain that \(\mathbf{w}_3+\mathbf{w}_6\) satisfies \[\begin{align} \label{eq-w3w6} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}({\boldsymbol{w}}_3+{\boldsymbol{w}}_6)=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{w}}_3+{\boldsymbol{w}}_6={\boldsymbol{\vartheta}}_3+{\boldsymbol{\vartheta}}_6, & \mathbf{x}\in\partial D,\\ {\boldsymbol{w}}_3+{\boldsymbol{w}}_6=\mathcal{O}(|\mathbf{x}|^{-1}),& |\mathbf{x}|\rightarrow\infty. \end{cases} \end{align}\tag{63}\] Multiplying the equation in 63 by \({\boldsymbol{w}}_3+{\boldsymbol{w}}_6\) and integrating by parts in \(B_R\setminus\overline{D}\), we have \[\begin{align} 2(\beta_{31}+\beta_{32})=-\int_{B_R\setminus\overline{D}}\left(\mathbb{C}e(\mathbf{w}_3+\mathbf{w}_6),e(\mathbf{w}_3+\mathbf{w}_6)\right)+\int_{\partial B_R}(\mathbf{w}_3+\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3+\mathbf{w}_6). \end{align}\] As \(|\mathbf{x}|\rightarrow\infty\), \[\label{R-condition} {\boldsymbol{w}}_3+{\boldsymbol{w}}_6=\mathcal{O}(|\mathbf{x}|^{-1})\quadand\quad \partial_{\boldsymbol{\nu}}({\boldsymbol{w}}_3+{\boldsymbol{w}}_6)=\mathcal{O}(|\mathbf{x}|^{-2}).\tag{64}\] This yields \[\begin{align} \label{int-w343w6} \int_{\partial B_R}(\mathbf{w}_3+\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3+\mathbf{w}_6)\rightarrow 0,\quadas~R\rightarrow\infty, \end{align}\tag{65}\] and thus, \[2(\beta_{31}+\beta_{32})=-\int_{D^e}\left(\mathbb{C}e(\mathbf{w}_3+\mathbf{w}_6),e(\mathbf{w}_3+\mathbf{w}_6)\right)<0.\]

In view of 27 and ?? again, we obtain \[\begin{align} \label{beta31-32} \beta_{31}-\beta_{32}=A_{33}-A_{36}=\int_{\partial D_1}\frac{\partial(\mathbf{w}_3-\mathbf{w}_6)}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\kappa}}_3=-\int_{\partial D_2}\frac{\partial(\mathbf{w}_3-\mathbf{w}_6)}{\partial{\boldsymbol{\nu}}}\Bigg|_{+}{\boldsymbol{\kappa}}_3. \end{align}\tag{66}\] Analogously to 63 , we have \[\begin{align} \label{eq-w3-w6} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}({\boldsymbol{w}}_3-{\boldsymbol{w}}_6)=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{w}}_3-{\boldsymbol{w}}_6={\boldsymbol{\vartheta}}_3-{\boldsymbol{\vartheta}}_6, & \mathbf{x}\in\partial D,\\ {\boldsymbol{w}}_3-{\boldsymbol{w}}_6={\boldsymbol{f}}_3-{\boldsymbol{f}}_6,& \mathbf{x}\in\partial B_R. \end{cases} \end{align}\tag{67}\] Multiplying the equation in 67 by \({\boldsymbol{w}}_3-{\boldsymbol{w}}_6\) and \({\boldsymbol{w}}_3+{\boldsymbol{w}}_6\), respectively, and integrating by parts in \(B_R\setminus\overline{D}\) yields \[\begin{align} \label{form-beta3132} 2(\beta_{31}-\beta_{32})=-\int_{B_R\setminus\overline{D}}\left(\mathbb{C}e(\mathbf{w}_3-\mathbf{w}_6),e(\mathbf{w}_3-\mathbf{w}_6)\right)+\int_{\partial B_R}(\mathbf{w}_3-\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3-\mathbf{w}_6) \end{align}\tag{68}\] and \[\label{w3-w6-w36} \int_{B_R\setminus\overline{D}}\left(\mathbb{C}e(\mathbf{w}_3-\mathbf{w}_6),e(\mathbf{w}_3+\mathbf{w}_6)\right)-\int_{\partial B_R}(\mathbf{w}_3+\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3-\mathbf{w}_6)=0,\tag{69}\] where we used 66 in the second equality. Similarly, multiplying the equation in 63 by \({\boldsymbol{w}}_3-{\boldsymbol{w}}_6\), integrating by parts in \(B_R\setminus\overline{D}\), and using 62 , we have \[\label{w3w6-w36} \int_{B_R\setminus\overline{D}}\left(\mathbb{C}e(\mathbf{w}_3+\mathbf{w}_6),e(\mathbf{w}_3-\mathbf{w}_6)\right)-\int_{\partial B_R}(\mathbf{w}_3-\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3+\mathbf{w}_6)=0.\tag{70}\] Applying 64 gives \[\begin{align} \label{w3-w6infty} \int_{\partial B_R}(\mathbf{w}_3-\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3+\mathbf{w}_6)\rightarrow 0\quadas~R\rightarrow\infty. \end{align}\tag{71}\] Since \((\mathbb{C}A,B)=(A,\mathbb{C}B)\) for any \(2\times2\) matrices \(A\) and \(B\) (see, for instance, [51]), it follows from 69 , 70 , and 71 that \[\begin{align} \int_{\partial B_R}(\mathbf{w}_3+\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3-\mathbf{w}_6)\rightarrow 0\quadas~R\rightarrow\infty. \end{align}\] This together with 71 and 65 implies that as \(R\rightarrow\infty\), the term \(\int_{\partial B_R}(\mathbf{w}_3-\mathbf{w}_6)\partial_{\boldsymbol{\nu}}(\mathbf{w}_3-\mathbf{w}_6)\) in 68 has the same asymptotics as \(4\int_{\partial B_R}\mathbf{w}_3\partial_{\boldsymbol{\nu}}\mathbf{w}_3\). Since the first term on the right-hand side of 68 is negative, it suffices to analyze \(\int_{\partial B_R}\mathbf{w}_3\partial_{\boldsymbol{\nu}}\mathbf{w}_3\) as \(R\rightarrow\infty\).

Recalling that \[\begin{align} \label{def-w3} {\boldsymbol{w}}_3:=\hat{\mathbf{S}}_D^\omega[\boldsymbol{\zeta}_3]=\int_{\partial D}(\mathbf{\Gamma}^{0,1}+ \mathbf{\Gamma}^{0,2})(\mathbf{x}-\mathbf{y})\boldsymbol{\zeta}_3(\mathbf{y})ds(\mathbf{y}). \end{align}\tag{72}\] It follows from 14 and 15 that \(\partial_{\boldsymbol{\nu}}\mathbf{\Gamma}^{0,1}=0\) and for \(|\mathbf{x}| \gg 1\) and \(\mathbf{y}\in \partial D\), we have \[\label{gamma02} \mathbf{\Gamma}^{0,2}(\mathbf{x}-\mathbf{y})=\frac{1}{4\pi}\Big(\frac{1}{\mu}+\frac{1}{\lambda +2\mu}\Big)\ln\left\vert\mathbf{x}\right\vert\mathbf{I}-\frac{1}{4\pi}\Big(\frac{1}{\mu}-\frac{1}{\lambda +2\mu}\Big)\frac{\mathbf{x}\otimes\mathbf{x}}{|\mathbf{x}|^2}+\mathcal{O}(|\mathbf{x}|^{-1}).\tag{73}\] By using 60 , we obtain on \(\partial B_R\), \[\begin{align} \partial_{\boldsymbol{\nu}}(\ln\left\vert\mathbf{x}\right\vert\mathbf{I}\boldsymbol{\zeta}_3(\mathbf{y}))=\frac{\mu}{R}\boldsymbol{\zeta}_3+\frac{\lambda+\mu}{R^3}(\boldsymbol{\zeta}_3\cdot\mathbf{x})\mathbf{x} \end{align}\] and \[\begin{align} \partial_{\boldsymbol{\nu}}\left(\frac{\mathbf{x}\otimes\mathbf{x}}{|\mathbf{x}|^2}\boldsymbol{\zeta}_3(\mathbf{y})\right)=\lambda \frac{(\boldsymbol{\zeta}_3\cdot\mathbf{x})\mathbf{x}}{R^3}+\mu\left(\frac{\boldsymbol{\zeta}_3}{R}-\frac{(\boldsymbol{\zeta}_3\cdot\mathbf{x})\mathbf{x}}{R^3}\right). \end{align}\] Thus, from 72 and 73 , we obtain \[\begin{align} \partial_{\boldsymbol{\nu}}\mathbf{w}_3&=a\int_{\partial D}\partial_{\boldsymbol{\nu}}(\ln\left\vert\mathbf{x}\right\vert\mathbf{I}\boldsymbol{\zeta}_3(\mathbf{y}))-b\int_{\partial D}\partial_{\boldsymbol{\nu}}\Big(\frac{\mathbf{x}\otimes\mathbf{x}}{|\mathbf{x}|^2}\boldsymbol{\zeta}_3(\mathbf{y})\Big)+o(1)\nonumber\\ &=\frac{\mu}{2\pi(\lambda+2\mu)R}I_{\boldsymbol{\zeta}_3}+\frac{1}{2\pi R^3}\Big(1+\frac{\lambda}{\lambda+2\mu}\Big)(I_{\boldsymbol{\zeta}_3}\cdot\mathbf{x})\mathbf{x}, \end{align}\] where \[a:=\frac{1}{4\pi}\Big(\frac{1}{\mu}+\frac{1}{\lambda +2\mu}\Big),\quad b:=\frac{1}{4\pi}\Big(\frac{1}{\mu}-\frac{1}{\lambda +2\mu}\Big),\] and \[I_{\boldsymbol{\zeta}_3}=(I_{\boldsymbol{\zeta}_3,1},I_{\boldsymbol{\zeta}_3,2}):=\int_{\partial D}\boldsymbol{\zeta}_3(\mathbf{y})ds(\mathbf{y}).\] Hence, we derive \[\begin{align} \label{int-w3-Dw3} &\int_{\partial B_R}\mathbf{w}_3\partial_{\boldsymbol{\nu}}\mathbf{w}_3\nonumber\\ &=\hat{\eta}_{\omega}\int_{\partial B_R} I_{\boldsymbol{\zeta}_3}\Big(\frac{\mu}{2\pi(\lambda+2\mu)R}I_{\boldsymbol{\zeta}_3}+\frac{1}{2\pi R^3}\big(1+\frac{\lambda}{\lambda+2\mu}\big)(I_{\boldsymbol{\zeta}_3}\cdot\mathbf{x})\mathbf{x}\Big)\nonumber\\ &\quad+a\int_{\partial B_R} (I_{\boldsymbol{\zeta}_3}\ln R)\Big(\frac{\mu}{2\pi(\lambda+2\mu)R}I_{\boldsymbol{\zeta}_3}+\frac{1}{2\pi R^3}\big(1+\frac{\lambda}{\lambda+2\mu}\big)(I_{\boldsymbol{\zeta}_3}\cdot\mathbf{x})\mathbf{x}\Big)\nonumber\\ &\quad-\frac{b}{R^2}\int_{\partial B_R} (I_{\boldsymbol{\zeta}_3}\mathbf{x}\otimes\mathbf{x})\Big(\frac{\mu}{2\pi(\lambda+2\mu)R}I_{\boldsymbol{\zeta}_3}+\frac{1}{2\pi R^3}\big(1+\frac{\lambda}{\lambda+2\mu}\big)(I_{\boldsymbol{\zeta}_3}\cdot\mathbf{x})\mathbf{x}\Big)\nonumber\\ &\quad+o(1)=:I+II+III+o(1). \end{align}\tag{74}\] Note that \[\int_{\partial B_R}(I_{\boldsymbol{\zeta}_3}\cdot\mathbf{x})^2=\int_0^{2\pi}(I_{\boldsymbol{\zeta}_3,1}R\cos\theta+I_{\boldsymbol{\zeta}_3,2}R\sin\theta)^2R\;d\theta=\pi R^3|I_{\boldsymbol{\zeta}_3}|^2,\] which yields \[\begin{align} \label{est-main-term} I&=\frac{\mu}{\lambda+2\mu}\hat{\eta}_{\omega}|I_{\boldsymbol{\zeta}_3}|^2+\frac{1}{2}\Big(1+\frac{\lambda}{\lambda+2\mu}\Big)\hat{\eta}_{\omega}|I_{\boldsymbol{\zeta}_3}|^2=\hat{\eta}_{\omega}|I_{\boldsymbol{\zeta}_3}|^2 \end{align}\tag{75}\] and \[\begin{align} II=\frac{1}{4\pi}\Big(\frac{1}{\mu}+\frac{1}{\lambda +2\mu}\Big)|I_{\boldsymbol{\zeta}_3}|^2\ln R. \end{align}\] Recalling that \(\beta_{3k}=\mathcal{O}(1)\) with \(k=1,2\), we note that \(\beta_{3k}\) should be independent of \(R\). Hence, the term \(II\) can be absorbed into the first term on the right-hand side of 68 . Direct calculations give \[\begin{align} \label{est-temrIII} III&=-\frac{b}{R^3}\Big(\frac{\mu}{2\pi(\lambda+2\mu)}+\frac{1}{2\pi}\big(1+\frac{\lambda}{\lambda+2\mu}\big)\Big)\int_{\partial B_R}(I_{\boldsymbol{\zeta}_3}\cdot\mathbf{x})^2\nonumber\\ &=-\frac{1}{8\pi}\Big(\frac{1}{\mu}-\frac{1}{\lambda +2\mu}\Big)\cdot\Big(1+\frac{\lambda+\mu}{\lambda+2\mu}\big)\Big)|I_{\boldsymbol{\zeta}_3}|^2. \end{align}\tag{76}\] Therefore, for any fixed \(|\omega|\ll 1\), as \(R\rightarrow\infty\), we derive from 68 , 74 , 75 , and 76 that \[\begin{align} \label{diff-beta3132} 2(\beta_{31}-\beta_{32})=\hat{\eta}_{\omega}|I_{\boldsymbol{\zeta}_3}|^2-\frac{1}{8\pi}\Big(\frac{1}{\mu}-\frac{1}{\lambda +2\mu}\Big)\cdot\Big(1+\frac{\lambda+\mu}{\lambda+2\mu}\big)\Big)|I_{\boldsymbol{\zeta}_3}|^2+\mathcal{O}(1), \end{align}\tag{77}\] where \(\mathcal{O}(1)\) is independent of \(\omega\).

We next analyze \(\hat{\eta}_{\omega}\). From 16 and 3 , it follows that \[\begin{align} \Re\hat{\eta}_{\omega}&=\frac{(\lambda+\mu)(16\pi\Re c_1-1)}{8\mu(\lambda+2\mu)\pi}+\frac{\Re\tilde{\eta}}{\mu}+\frac{\ln|\omega|-\ln{c_s}}{4\mu\pi}-\frac{\ln|\omega|-\ln{c_p}}{4(\lambda+2\mu)\pi}\\ &=\frac{(\lambda+\mu)(16\pi\Re c_1-1)}{8\mu(\lambda+2\mu)\pi}+\frac{\Re\tilde{\eta}}{\mu}+\frac{(\lambda+\mu)\ln|\omega|}{4\mu(\lambda+2\mu)\pi}-\frac{\ln{c_s}}{4\mu\pi}+\frac{\ln{c_p}}{4(\lambda+2\mu)\pi}, \end{align}\] and \[\begin{align} \Im\hat{\eta}_{\omega}&=\frac{16\pi(\lambda+\mu)\Im c_1}{8\mu(\lambda+2\mu)\pi}+\frac{\Im\tilde{\eta}}{\mu}+\frac{(\lambda+\mu)\theta}{4\mu(\lambda+2\mu)\pi}. \end{align}\] Here, we used \(\omega=|\omega|e^{\mathrm{i}\theta}\) with \(\theta=\arg\omega\). It is well known established that resonant frequencies are located in the lower half of the complex plane [61]. Consequently, applying the conditions \(|\omega|\ll 1\) and \(\Re\omega> 0\) from Definition 1, we deduce that \(\theta\in(-\pi/2,0)\). Therefore, \[\Re\hat{\eta}_{\omega}<0,\quad \Im\hat{\eta}_{\omega}<0.\] Coming back to 77 , we derive \(\Re(\beta_{31}-\beta_{32})<0\) and \(\Im(\beta_{31}-\beta_{32})<0\) when \(|\omega|\ll 1\). This completes the proof of Proposition 2. ◻

Proof. By the classical gradient estimates for elliptic systems (see, for instance, [60]), we have \[\label{est-u11-out} \|\nabla \mathbf{u}_{i}\|_{L^\infty(\mathbb{R}^2\setminus\overline{B}_R)}\leq C,\quad i=1,\dots,6,\tag{83}\] where \(C>0\) is a constant independently of \(\varepsilon\). Thus, we shall prove the estimates of \(|\nabla \mathbf{u}_{i}(\mathbf{x})|\) for \(\mathbf{x}\in B_R\setminus\overline{D}=:\Omega\).

Step 1: Estimate of \(\nabla{\mathbf{u}}_i\) with \(i=1,5,6\). We provide the details for the case \(i=1\) since the cases \(i=5,6\) are analogous and thus omitted. Decompose \({\mathbf{u}}_1\) as \[\label{dec-u1} {\mathbf{u}}_1={\mathbf{u}}_{1,1}+{\mathbf{u}}_{1,2},\tag{84}\] where \[\begin{align} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}{\boldsymbol{u}}_{1,1}+\rho\omega_1^2{\boldsymbol{u}}_{1,1}=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{u}}_{1,1}={\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{1,1}, & \mathbf{x}\in\partial D_1,\\ {\boldsymbol{u}}_{1,1}=-{\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{1,2}, & \mathbf{x}\in\partial D_2,\\ {\boldsymbol{u}}_{1,1}=0,& \mathbf{x}\in\partial B_R, \end{cases} \end{align}\] and \[\begin{align} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}{\boldsymbol{u}}_{1,2}+\rho\omega_1^2{\boldsymbol{u}}_{1,2}=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{u}}_{1,2}=0, & \mathbf{x}\in\partial D_1,\\ {\boldsymbol{u}}_{1,2}=0, & \mathbf{x}\in\partial D_2,\\ {\boldsymbol{u}}_{1,2}={\boldsymbol{g}}_1,& \mathbf{x}\in\partial B_R. \end{cases} \end{align}\] Applying a similar argument in the proof of Lemma 14, we obtain \[\label{est-Du12} |\nabla\mathbf{u}_{1,2}(\mathbf{x})|\leq C\big(\|\mathbf{u}_{1,2}\|_{L^{2}(D^e)}+\|{\boldsymbol{g}}_1\|_{L^{\infty}(\partial B_R)}\big),\quad \mathbf{x}\in D^{e}.\tag{85}\] To estimate \(\nabla{\boldsymbol{u}}_{1,1}\), we construct an auxiliary function \(\bar{\boldsymbol{u}}_{1,1}\in C^{2,\alpha}(\mathbb{R}^2)\), such that \(\bar{\boldsymbol{u}}_{1,1}={\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{1,1}\) on \(\partial D_1\), \(\bar{\boldsymbol{u}}_{1,1}=-{\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{1,2}\) on \(\partial D_2\), \(\bar{\boldsymbol{u}}_{1,1}=0\) on \(\partial B_R\), \[\label{def-baru11} \bar{\boldsymbol{u}}_{1,1}=p({\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{1,1})+(1-p)(-{\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{1,2})\quadin\;\Omega_{2R_0},\tag{86}\] and \(\|\bar {\boldsymbol{u}}_{1,1}\|_{C^{2,\alpha}(D^e\setminus\Omega_{R_0})}\leq\,C\), where \(p\) is defined in 56 . Using some calculations, we have \[\label{est-baru11} |\partial_{{\boldsymbol{x}}_1}\bar {\boldsymbol{u}}_{1,1}^{(2)}|\leq \frac{C|{\boldsymbol{x}}_1|}{\delta({\boldsymbol{x}}_1)}\quadand\quad \partial_{{\boldsymbol{x}}_2}\bar {\boldsymbol{u}}_{1,1}^{(2)}=\frac{2+o(1)}{\delta({\boldsymbol{x}}_1)}.\tag{87}\] Denote \[\bar {\boldsymbol{w}}_{1,1}:={\boldsymbol{u}}_{1,1}-\bar {\boldsymbol{u}}_{1,1},\] then we have \[\begin{align} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}\bar{\boldsymbol{w}}_{1,1}+\rho\omega_1^2\bar{\boldsymbol{w}}_{1,1}=-\mathcal{L}_{ {\lambda}, {\mu}}\bar{\boldsymbol{u}}_{1,1}-\rho\omega_1^2\bar{\boldsymbol{u}}_{1,1}, & \mathbf{x}\in D^{e},\\ {\boldsymbol{w}}_{1,1}=0, & \mathbf{x}\in\partial D,\\ {\boldsymbol{w}}_{1,1}=0,& \mathbf{x}\in\partial B_R. \end{cases} \end{align}\] By applying a similar argument in the proof of [54] (see also, Lemma 13), we have \[\label{est-u11-} \nabla{\boldsymbol{u}}_{1,1}=\nabla\bar{\boldsymbol{u}}_{1,1}+\mathcal{O}(1).\tag{88}\] Thus, combining with 86 and 87 , we obtain \[\label{est-u11} |\nabla{\boldsymbol{u}}_{1,1}(\mathbf{x})|\leq \frac{C}{\varepsilon+\kappa|{\boldsymbol{x}}_1|^2}+C\big(\|{\boldsymbol{f}}_{1,1}\|_{C^2(\partial D_1)}+\|{\boldsymbol{f}}_{1,2}\|_{C^2(\partial D_2)}\big),\quad\mathbf{x}\in\Omega_{R_0},\tag{89}\] and \[|\nabla{\boldsymbol{u}}_{1,1}(0,{\boldsymbol{x}}_2)|\geq \frac{1}{C\varepsilon}.\] Combining 83 , 84 , 85 , and 89 , we conclude \[\begin{align} |\nabla{\boldsymbol{u}}_{1}(\mathbf{x})|&\leq \frac{C}{\varepsilon+\kappa|{\boldsymbol{x}}_1|^2}+C\big(\|\mathbf{u}_{1}\|_{L^2(D^e)}+\|{\boldsymbol{g}}_{1}\|_{L^{\infty}(\partial B_R)}+\|{\boldsymbol{f}}_{1,1}\|_{C^2(\partial D_1)}\\ &\quad+\|{\boldsymbol{f}}_{1,2}\|_{C^2(\partial D_2)}\big),\quad\mathbf{x}\in D^e. \end{align}\]

Step 2: Estimate of \(\nabla{\mathbf{u}}_2\). As in 84 , we decompose \({\mathbf{u}}_2\) into \[{\mathbf{u}}_2={\mathbf{u}}_{2,1}+{\mathbf{u}}_{2,2},\] where \[\begin{align} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}{\boldsymbol{u}}_{2,1}=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{u}}_{2,1}={\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{2,1}, & \mathbf{x}\in\partial D_1,\\ {\boldsymbol{u}}_{2,1}={\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{2,2}, & \mathbf{x}\in\partial D_2,\\ {\boldsymbol{u}}_{2,1}=0,& \mathbf{x}\in\partial B_R, \end{cases} \end{align}\] and \[\begin{align} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}{\boldsymbol{u}}_{2,2}=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{u}}_{2,2}=0, & \mathbf{x}\in\partial D_1,\\ {\boldsymbol{u}}_{2,2}=0, & \mathbf{x}\in\partial D_2,\\ {\boldsymbol{u}}_{2,2}={\boldsymbol{g}}_2,& \mathbf{x}\in\partial B_R. \end{cases} \end{align}\] Then by using similar arguments that led to Lemma 14, we have \[|\nabla{\boldsymbol{u}}_{2,2}(\mathbf{x})|\leq C\big(\|\mathbf{u}_{2,2}\|_{L^2(D^e)}+\|{\boldsymbol{g}}_{2}\|_{L^{\infty}(\partial B_R)}\big),\quad\mathbf{x}\in D^e.\] In order to prove the estimate of \(\nabla{\boldsymbol{u}}_{2,1}\), we use the auxiliary function \(p\) in 56 to construct \(\bar{\boldsymbol{u}}_{2,1}\in C^{2,\alpha}(\mathbb{R}^2)\), such that \(\bar{\boldsymbol{u}}_{2,1}={\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{2,1}\) on \(\partial D_1\), \(\bar{\boldsymbol{u}}_{2,1}={\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{2,2}\) on \(\partial D_2\), \(\bar{\boldsymbol{u}}_{2,1}=0\) on \(\partial B_R\), \[\bar{\boldsymbol{u}}_{2,1}=p({\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{2,1})+(1-p)({\boldsymbol{\kappa}}_2+{\boldsymbol{f}}_{2,2})\quadin\;\Omega_{2R_0},\] and \(\|\bar {\boldsymbol{u}}_{2,1}\|_{C^{2,\alpha}(D^e\setminus\Omega_{R_0})}\leq\,C\). Then similar to 88 , we have \[\nabla{\boldsymbol{u}}_{2,1}=\nabla\bar{\boldsymbol{u}}_{2,1}+\mathcal{O}(1),\] which gives \[\begin{align} |\nabla{\boldsymbol{u}}_{2,1}(\mathbf{x})|&\leq \frac{C|{\boldsymbol{f}}_{2,1}({\boldsymbol{x}}_1,\varepsilon/2+\mathcal{H}_1({\boldsymbol{x}}_1))-{\boldsymbol{f}}_{2,2}({\boldsymbol{x}}_1,-\varepsilon/2+\mathcal{H}_2({\boldsymbol{x}}_1))|}{\varepsilon+\kappa|{\boldsymbol{x}}_1|^2}\\ &\quad+C\big(\|{\boldsymbol{f}}_{2,1}\|_{C^2(\partial D_1)}+\|{\boldsymbol{f}}_{2,2}\|_{C^2(\partial D_2)}\big). \end{align}\] Therefore, we obtain \[\begin{align} |\nabla{\boldsymbol{u}}_{2}(\mathbf{x})|&\leq \frac{C|{\boldsymbol{f}}_{2,1}({\boldsymbol{x}}_1,\varepsilon/2+\mathcal{H}_1({\boldsymbol{x}}_1))-{\boldsymbol{f}}_{2,2}({\boldsymbol{x}}_1,-\varepsilon/2+\mathcal{H}_2({\boldsymbol{x}}_1))|}{\varepsilon+\kappa|{\boldsymbol{x}}_1|^2}\\ &\quad+C\big(\|\mathbf{u}_{2}\|_{L^2(D^e)}+\|{\boldsymbol{g}}_{2}\|_{L^{\infty}(\partial B_R)}+\|{\boldsymbol{f}}_{2,1}\|_{C^2(\partial D_1)}+\|{\boldsymbol{f}}_{2,2}\|_{C^2(\partial D_2)}\big). \end{align}\]

Step 3: Estimates of \(\nabla{\mathbf{u}}_i\) with \(i=3,4\). By following similar arguments in Steps 1 and 2, we decompose \({\mathbf{u}}_i\) into \[{\mathbf{u}}_i={\mathbf{u}}_{i,1}+{\mathbf{u}}_{i,2},\] where \[\begin{align} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}{\boldsymbol{u}}_{i,1}=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{u}}_{i,1}=e_{i,1}{\boldsymbol{\kappa}}_1+e_{i,3}{\boldsymbol{\kappa}}_3+{\boldsymbol{f}}_{i,1}, & \mathbf{x}\in\partial D_1,\\ {\boldsymbol{u}}_{i,1}=e_{i,1}{\boldsymbol{\kappa}}_1+e_{i,6}{\boldsymbol{\kappa}}_3+{\boldsymbol{f}}_{i,2}, & \mathbf{x}\in\partial D_2,\\ {\boldsymbol{u}}_{i,1}=0,& \mathbf{x}\in\partial B_R, \end{cases} \end{align}\] and \[\begin{align} \begin{cases} \mathcal{L}_{ {\lambda}, {\mu}}{\boldsymbol{u}}_{i,2}=0, & \mathbf{x}\in D^{e},\\ {\boldsymbol{u}}_{i,2}=0, & \mathbf{x}\in\partial D_1,\\ {\boldsymbol{u}}_{i,2}=0, & \mathbf{x}\in\partial D_2,\\ {\boldsymbol{u}}_{i,2}={\boldsymbol{g}}_i,& \mathbf{x}\in\partial B_R,\quad i=3,4. \end{cases} \end{align}\] Then we have \[|\nabla{\boldsymbol{u}}_{i,2}(\mathbf{x})|\leq C\big(\|\mathbf{u}_{i,2}\|_{L^2(D^e)}+\|{\boldsymbol{g}}_{i}\|_{L^{\infty}(\partial B_R)}\big),\quad\mathbf{x}\in D^e,\quad i=3,4.\] We next construct \(\bar{\boldsymbol{u}}_{i,1}\in C^{2,\alpha}(\mathbb{R}^2)\), such that \(\bar{\boldsymbol{u}}_{i,1}=e_{i,1}{\boldsymbol{\kappa}}_1+e_{i,3}{\boldsymbol{\kappa}}_3+{\boldsymbol{f}}_{i,1}\) on \(\partial D_1\), \(\bar{\boldsymbol{u}}_{i,1}=e_{i,1}{\boldsymbol{\kappa}}_1+e_{i,6}{\boldsymbol{\kappa}}_3+{\boldsymbol{f}}_{i,2}\) on \(\partial D_2\), \(\bar{\boldsymbol{u}}_{i,1}=0\) on \(\partial B_R\), \[\bar{\boldsymbol{u}}_{i,1}=p(e_{i,1}{\boldsymbol{\kappa}}_1+e_{i,3}{\boldsymbol{\kappa}}_3+{\boldsymbol{f}}_{i,1})+(1-p)(e_{i,1}{\boldsymbol{\kappa}}_1+e_{i,6}{\boldsymbol{\kappa}}_3+{\boldsymbol{f}}_{i,2})\quadin\;\Omega_{2R_0},\] and \(\|\bar {\boldsymbol{u}}_{i,1}\|_{C^{2,\alpha}(D^e\setminus\Omega_{R_0})}\leq\,C\). Direct computations give \[\partial_{{\boldsymbol{x}}_1}\bar {\boldsymbol{u}}_{i,1}^{(1)}=\frac{|{\boldsymbol{x}}_1|}{\delta({\boldsymbol{x}}_1)}o(1)+|x_1|,\quad\partial_{{\boldsymbol{x}}_2}\bar {\boldsymbol{u}}_{i,1}^{(1)}=\frac{1}{\delta({\boldsymbol{x}}_1)}o(1)+\mathcal{O}(1),\] and \[\partial_{{\boldsymbol{x}}_1}\bar {\boldsymbol{u}}_{i,1}^{(2)}=\frac{|{\boldsymbol{x}}_1|}{\delta({\boldsymbol{x}}_1)}o(1)+\mathcal{O}(1),\quad \partial_{{\boldsymbol{x}}_2}\bar {\boldsymbol{u}}_{i,1}^{(2)}=\frac{2\mathbf{x}_1+o(1)}{\delta({\boldsymbol{x}}_1)}.\] Similar to 88 , we have \[\nabla{\boldsymbol{u}}_{i,1}=\nabla\bar{\boldsymbol{u}}_{i,1}+\mathcal{O}(1),\quad i=3,4.\] Thus, we obtain \[\begin{align} |\nabla{\boldsymbol{u}}_{i}(\mathbf{x})|&\leq C\left(\frac{|{\boldsymbol{x}}_1|}{\varepsilon+\kappa|{\boldsymbol{x}}_1|^2}+1\right)+C\big(\|\mathbf{u}_{i}\|_{L^2(D^e)}+\|{\boldsymbol{g}}_{i}\|_{L^{\infty}(\partial B_R)}+\|{\boldsymbol{f}}_{i,1}\|_{C^2(\partial D_1)}\\ &\quad+\|{\boldsymbol{f}}_{i,2}\|_{C^2(\partial D_2)}\big),\quad i=3,4,\quad\mathbf{x}\in D^e. \end{align}\] When \(|{\boldsymbol{x}}_1|\leq\varepsilon\), we get that \(|\nabla{\boldsymbol{u}}_{i}(\mathbf{x})|\) with \(i=3,4\) are bounded. When \(|{\boldsymbol{x}}_1|>\varepsilon\), \(|\nabla{\boldsymbol{u}}_{i}(\mathbf{x})|\) can blow up and \[|\nabla{\boldsymbol{u}}_{i}({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)|\geq|\partial_{{\boldsymbol{x}}_2}{\boldsymbol{u}}_{i}^{(2)}({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)| \geq\frac{1}{C\sqrt{\varepsilon}},\quad |{\boldsymbol{x}}_1|=\sqrt{\kappa^{-1}\varepsilon},\quad i=3,4.\] We complete the proof. ◻