\section{Introduction}

During last three decades or so, there has been significant progress in quantitative analysis of the field concentration (field enhancement or stress) in the narrow region between two inclusions for the homogeneous and inhomogeneous transmission problems of conductivity equations (see, for example, \cite{AKLLL, AKL, BLY, DL, LV, Yun}). In this paper we revisit the problem by exploiting the spectral nature residing behind the phenomenon of the field concentration when inclusions are circular. Through the spectral approach we not only gain new insights on the existing results but also obtain significant new results. We show that when $(k_1-1)(k_2-1) <0$, where $k_1, \ k_2$ are conductivities of inclusions and $1$ is conductivity of the background, have different signs (typically, $k_1=0$ and $k_2=\infty$), then the gradient of the solution is bounded regardless of the distance between inclusions, but the second and higher derivatives may blow up as the distance tends to zero. We show that this result holds for both homogeneous and inhomogeneous problems.

We now state the main results of this paper in a precise manner and then compare them with existing results. Let $D_1$ and $D_2$ be disks in $\mathbb{R}^2$ which are separated by the small distance denoted by $\epsilon$. Let $D:=D_1 \cup D_2$. Let $k_j$ be the conductivity of $D_j$ for $j=1,2$, while that of $\mathbb{R}^2 \setminus D$ is assumed to be $1$. So the conductivity distribution is given by
\begin{equation}
\sigma = \chi_{\mathbb{R}^2 \setminus D} + k_1 \chi_{D_1}  + k_2 \chi_{D_2},
\end{equation}
where $\chi$ denotes the characteristic function on the respective set. Assuming that $0<k_j \neq 1 <\infty$ ($j=1,2$), we consider the inhomogeneous transmission problem: for a given function $f$
\begin{equation}\label{NCE}
\begin{cases}
	\nabla \cdot \sigma \nabla u = f \quad &\text{in } \mathbb{R}^2, \\
	u(x) = c \ln |x| + O(|x|^{-1}) \quad &\text{as } |x| \to \infty
\end{cases}
\end{equation}
for some constant $c$. We also consider the homogeneous transmission problem:
\begin{equation}\label{HCE}
\begin{cases}
\nabla \cdot \sigma \nabla u = 0 \quad &\text{in } \mathbb{R}^2, \\
u(x) - H(x) = O(|x|^{-1}) \quad &\text{as } |x| \to \infty,
\end{cases}
\end{equation}
where $H$ is a given function harmonic in $\mathbb{R}^2$.

Note that there are three parameters involved in the problem \eqnref{NCE}: conductivities $k_1$, $k_2$, and the distance $\epsilon$ between two inclusions. The problem is to derive estimates for derivatives of $u$ in terms of these parameters when $\epsilon$ tends to $0$.

We assume that the inhomogeneity $f$ is given by $f=\nabla \cdot g$ for some $g$. It is assumed that $g$ is compactly supported in $\mathbb{R}^2$, which is just for the sake of simplicity and can be replaced with some integrability condition. We also impose some regularity condition on $g$: $g$ is piecewise $C^{n, \alpha}$ for some non-negative integer $n$ and $0<\alpha<1$. The function $g$ being piecewise $C^{n, \alpha}$ in this paper means that $g$ is $C^{n,\alpha}$ on $\overline{D_1}$, $\overline{D_2}$ and $\mathbb{R}^2 \setminus D$ separately. For piecewise $C^{n, \alpha}$ functions $g$, the norm is defined by
\begin{equation}\label{Cknorm}
\| g \|_{n,\alpha}:= \|g\|_{C^{n,\alpha}(\overline{D_1})} +  \|g\|_{C^{n,\alpha}(\overline{D_2})}+ \|g\|_{C^{n,\alpha}(\mathbb{R}^2 \setminus D)}.
\end{equation}
When $\alpha=0$, we denote it by $\| g \|_{n,0}$. We also use the following norm:
\begin{equation}\label{Ckstarnorm}
\| g \|_{n,\alpha}^*:= \frac{1}{k_1} \|g\|_{C^{n,\alpha}(\overline{D_1})} +  \frac{1}{k_2} \|g\|_{C^{n,\alpha}(\overline{D_2})}+ \|g\|_{C^{n,\alpha}(\mathbb{R}^2 \setminus D)}.
\end{equation}
We also denote by $\| u \|_{n,\Omega}$ the piece-wise $C^n$ norm on $\Omega$ when $\Omega$ is a bounded set containing $\overline{D_1 \cup D_2}$.

When $(k_1-1)(k_2-1) >0$, we obtain the following theorems for the inhomogeneous and homogeneous transmission problems. Here and throughout this paper, we put
\begin{equation}\label{notation2}
r_*:= \sqrt{\frac{2(r_1+r_2)}{r_1 r_2}}
\end{equation}
for ease of notation ($r_j$ is the radius of $D_j$, $j=1,2$). We also put
\begin{equation}\label{Glj}
\lambda_j= \frac{k_j+1}{2(k_j-1)}, \quad j=1,2.
\end{equation}

\begin{theorem}\label{thm:main1}
Suppose $(k_1-1)(k_2-1) >0$ and $f=\nabla \cdot g$ for some piecewise $C^{n-1, \alpha}$ function $g$ with the compact support ($n$ is a positive integer and $0<\alpha<1$). There is a constant $C>0$ independent of $k_1$, $k_2$, $\epsilon$ and $g$ such that the solution $u$ to \eqnref{NCE} satisfies
\begin{equation}\label{mainest1}
\| u \|_{n,0} \le C \| g \|_{n-1,\alpha}^* \left( 4\lambda_1\lambda_2-1 + r_* \sqrt{\epsilon} \right)^{-n}.
\end{equation}
This estimate is optimal in the sense that there is $f$ such that the reverse inequality (with a different constant $C$) holds when $n=1$.
\end{theorem}

\begin{theorem}\label{thm:main4}
Let $\Omega$ be a bounded set containing $\overline{D_1 \cup D_2}$. Let $u$ be the solution to \eqnref{HCE}. If $(k_1-1)(k_2-1) >0$, then there is a constant $C>0$ independent of $k_1$, $k_2$, $\epsilon$ and the function $H$ such that
\begin{equation}\label{mainest4-1}
\| u \|_{n,\Omega} \le C \| H \|_{C^n (\Omega)} \left( 4\lambda_1\lambda_2-1 + r_* \sqrt{\epsilon} \right)^{-n} .
\end{equation}
This estimate is optimal in the sense that there is a harmonic function $H$ such that the reverse inequality (with a different constant $C$) holds for the case $n=1$.
\end{theorem}

Since
$$
4\lambda_1\lambda_2-1 = \frac{2(k_1+k_2)}{(k_1-1)(k_2-1)},
$$
the estimate \eqnref{mainest1} shows that if either $k_1$ or $k_2$ is finite (away from $0$ and $\infty$), then $\| u \|_{n,0}$ is bounded regardless of the distance $\epsilon$, while if both $k_1$ and $k_2$ tend to $\infty$, then the right-hand side of \eqnref{mainest1} is of order ${\epsilon}^{-n/2}$. As shown by an example in section \ref{sec:opti}, $\nabla u$ may actually blow up at the order of ${\epsilon}^{-1/2}$ . If $k_1$ and $k_2$ tend to $0$, then the right-hand side of \eqnref{mainest1} is also of order ${\epsilon}^{-n/2}$ provided that $\| g \|_{n-1,\alpha}^*$ is bounded, in particular, if there is no source in $D_1 \cup D_2$, namely, $g=0$ in $D_1 \cup D_2$. The estimate \eqnref{mainest4-1} yields the same findings.

The estimates \eqnref{mainest1} and \eqnref{mainest4-1} are not new; Here, we derive them using a different method: the spectrum of the Neumann-Poincar\'e operator related to the problems. The estimate \eqnref{mainest1} was obtained in \cite{DL}.
The estimate \eqnref{mainest4-1} has a long history: The estimate for the gradient, namely, for $n=1$, is obtained in \cite{AKLLL, AKL}; that for higher $n$ in \cite{DL}. The results for circular inclusions have been extended to those of more general shape with the conductivity $\infty$ \cite{BLY, KY19, Lekner-PRSA-12, LY-CPDE-09, Yun}, the insulating problem in three dimensions or higher \cite{BLY10, LY, Wein, Yun16}, and to other equations; $p-$Laplacian \cite{GN-MMS-12}, Lam\'e system \cite{BLL-ARMA-15, BLL-AM-17, KY}, and Stokes system \cite{AKKY}. The estimate \eqnref{mainest4-1} says in particular that $\nabla u$ is bounded regardless of the distance between two inclusions if $k_1$ or $k_2$ is finite. This fact is known to be true in a more general setting \cite{BV, LN, LV}. Asymptotic characterizations of the gradient blow-up when $k_1=k_2=\infty$ (or $k_1=k_2=0$) are obtained in \cite{ACKLY13, KLY, KLY15, LY-JMAA-15}.

Theorem \ref{thm:main1} and \ref{thm:main4} (and theorems to follow) are proved by exploiting the spectral nature of the problem \eqnref{NCE}. The problem \eqnref{NCE} can be reduced to the following problem (see the beginning of Section \ref{sec:NP}):
\begin{equation}\label{mainproblem}
\begin{cases}
\Delta v = 0  \quad&\mbox{in } D \cup (\mathbb{R}^2 \setminus \overline{D}), \\
\displaystyle v|_+ - v|_- = 0  \quad&\mbox{on }\partial D,  \\
\displaystyle \partial_\nu v|_+ - k_j \partial_\nu v|_- = (k_j-1)\eta_j  \quad&\mbox{on }\partial D_j, \ j=1,2, \\
\displaystyle v (x) = c \ln |x| + O( |x|^{-1})~&\mbox{as }|x|\rightarrow\infty
\end{cases}
\end{equation}
for some constant $c$. Here and throughout this paper, $\partial_\nu$ denotes the outward normal derivative on $\partial D_j$ and the subscripts $\pm$ denote the limits from outside and inside of $D_j$ respectively. . For the problem \eqnref{NCE}, $\eta_j$ appearing in \eqnref{mainproblem} is given by $\eta_j=\partial_\nu F|_{\partial D_j}$ ($j=1,2$), where $F$ is the (weighted) Newtonian potential of $f$, namely,
\begin{equation}\label{Fdef}
F(x)=  \frac{1}{2\pi} \int_{\mathbb{R}^2 \setminus D} \ln|x-y| f(y) \, dy + \sum_{j=1}^2 \frac{1}{2 \pi k_j} \int_{D_j} \ln|x-y| f(y) \, dy
\end{equation}
for $x \in \mathbb{R}^2$. For the problem \eqnref{HCE}, $\eta_j=\partial_\nu H|_{\partial D_j}$.

If we represent the solution to \eqnref{mainproblem} in terms of the single layer potential, then the interface condition on $\partial D_j$ for $j=1,2$ (the third line in \eqnref{mainproblem}) is reduced to an integral equation for the Neumann-Poincar\'e (abbreviated by NP) operator on two circles $\partial D_1$ and $\partial D_2$ (see section \ref{sec:NP}). By solving the integral equation we see explicitly that the gradient and higher derivatives of the solution to \eqnref{mainproblem} is quantitatively estimated in terms of the quantity $4 \lambda_1 \lambda_2-\rho$, where $\rho$ is a constant determined by the radii $r_1$ and $r_2$ of inclusions $D_1$ and $D_2$ and the distance $\epsilon$ between them such that $\pm \frac{1}{2} \rho^n$ ($n=1, 2, \ldots$) are eigenvalues of the NP operator (see section \ref{sec:twocircles} and \ref{sec:specrep}). Since $\rho \sim 1- r_* \sqrt{\epsilon}$ (see \eqnref{Grasymp}), we have estimates \eqnref{mainest1} and \eqnref{mainest4-1}. Here and throughout this paper $A \sim B$ for some positive quantities $A$ and $B$ depending on some parameters means that there are positive constants $C_1$ and $C_2$ independent of those parameters such that $C_1 \le A/B \le C_2$. Notation $\lesssim$ and $\gtrsim$ are defined analogously.

If $(k_1-1)(k_2-1) <0$, then $4 \lambda_1 \lambda_2 <0$. Thus the right-hand sides of \eqnref{mainest1} and \eqnref{mainest4-1} are bounded and are not optimal bounds in this case. We obtain the following theorem for this case.

\begin{theorem}\label{thm:main3}
Suppose $(k_1-1)(k_2-1) <0$ and $f=\nabla \cdot g$ for some piecewise $C^{n, \alpha}$ function $g$ with the compact support ($n$ is a positive integer and $0<\alpha<1$). There is a constant $C>0$ independent of $k_1$, $k_2$, $\epsilon$ and $g$ such that the solution $u$ to \eqnref{NCE} satisfies
\begin{equation}\label{mainest3-1}
\| u \|_{n,0} \le C \| g \|_{n,\alpha}^* \left( 4|\lambda_1 \lambda_2|-1 + r_* \sqrt{\epsilon} \right)^{-n+1}.
\end{equation}
This estimate is optimal in the sense that there is $f$ such that the reverse inequality (with a different constant $C$) holds for $n=2$.
\end{theorem}

The estimate \eqnref{mainest3-1} shows that if $(k_1-1)(k_2-1) <0$, then $\nabla u$ is bounded regardless of the $k_1$, $k_2$ and $\epsilon$, provided that $\| g \|_{n,\alpha}^*$ is bounded. But, the $n$th ($n \ge 2$) order derivative may blow up at the rate of $\epsilon^{-(n-1)/2}$ if, for example, $k_1=0$ and $k_2=\infty$. The second derivative of $u$ actually blows up at the rate of $\epsilon^{-1/2}$ in some case as we show in section \ref{sec:opti}.

We obtain the following theorem for the homogeneous problem similar to Theorem \ref{thm:main3}.

\begin{theorem}\label{thm:main5}
Let $\Omega$ be a bounded set containing $\overline{D_1 \cup D_2}$. Let $u$ be the solution to \eqnref{HCE}. If $(k_1-1)(k_2-1) <0$, then there is a constant $C>0$ independent of $k_1$, $k_2$, $\epsilon$ and the function $H$ such that
\begin{equation}\label{mainest4-2}
\| u \|_{n,\Omega} \le C \| H \|_{C^{n+1} (\Omega)} \left( 4|\lambda_1 \lambda_2|-1 + r_* \sqrt{\epsilon} \right)^{-n+1}.
\end{equation}
This estimate is optimal in the sense that there is a harmonic function $H$ such that the reverse inequality (with a different constant $C$) holds for $n=2$.
\end{theorem}

This paper is organized as follows. In section \ref{sec:NP}, we review the spectral theory of the NP operator defined on the boundary of two disjoint inclusions of general shape and prove uniqueness and existence of the solution to \eqnref{mainproblem}. Section \ref{sec:twocircles} is to show M\"obius invariance of the NP spectra and derive the full spectrum on two circles. In section \ref{sec:specrep} we derive representation formulas for the solutions to \eqnref{NCE} and \eqnref{HCE} in terms of NP spectrum. Section \ref{sec:pfs} is to prove Theorem \ref{thm:main1}-\ref{thm:main5} and section \ref{sec:opti} is to show optimality of estimates in those theorems.

\section{The NP operator and well-posedness of the problem}\label{sec:NP}

Let us begin this section by showing how the problem \eqnref{NCE} is reduced to \eqnref{mainproblem}. The problem \eqnref{NCE} can be expressed as
\begin{equation}\label{200}
\begin{cases}
\Delta u = f  \quad&\mbox{in } \mathbb{R}^2 \setminus \overline{D}, \\
\Delta u = k_j^{-1} f  \quad&\mbox{in } D_j, \ j=1,2 , \\
\displaystyle u|_+ - u|_- = 0  \quad&\mbox{on }\partial D \\
\displaystyle \partial_\nu u|_+ - k_j \partial_\nu u|_- = 0  \quad&\mbox{on }\partial D_j, \ j=1,2 ,\\
\displaystyle u (x) = c \ln|x| + O( |x|^{-1})~&\mbox{as }|x|\rightarrow\infty.
\end{cases}
\end{equation}
Let $v:=u-F$ with $F$ given by \eqnref{Fdef}. Since $\Delta F= f$ in $\mathbb{R}^2 \setminus \overline{D}$ and $\Delta F = k_j^{-1} f$ in $D_j$, $v$ is the solution to \eqnref{mainproblem}. By putting $v:u-H$, we see that \eqnref{HCE} is reduced to \eqnref{mainproblem}.

In this section we show that the problem \eqnref{mainproblem} is reduced to an integral equation of the NP operator. We also prove uniqueness and existence of the solution to \eqnref{mainproblem} for $\eta_j \in H^{-1/2} (\partial D_j)$ for $j=1,2$, and of \eqnref{NCE} as a consequence. Here and afterwards, $H^{-1/2} (\partial D_j)$ denotes the Sobolev space of order $-1/2$ on $\partial D_j$, and $H^{-1/2}_0 (\partial D_j)$ is the subspace of $H^{-1/2} (\partial D_j)$ whose element $\eta_j$ satisfies
\begin{equation}
\int_{\partial D_j} \eta_j d\sigma =0.
\end{equation}
We use $\mathcal{H}^{-1/2}=\mathcal{H}^{-1/2}(\partial D)$ and $\mathcal{H}^{-1/2}_0=\mathcal{H}_0^{-1/2}(\partial D)$ respectively to denote $H^{-1/2} (\partial D_1) \times H^{-1/2} (\partial D_2)$ and $H^{-1/2}_0 (\partial D_1) \times H^{-1/2}_0 (\partial D_2)$ for ease of notation.

\subsection{Neumann-Poincar\'e operator}

Let $\Omega$ be a simply connected bounded domain in $\mathbb{R}^2$ with the Lipschitz continuous boundary. For a function $\mu \in H^{-1/2}(\partial \Omega)$, the single layer potential $\mathcal{S}_{\partial \Omega} [\mu]$ is defined by
\begin{equation}
\mathcal{S}_{\partial \Omega} [\mu] (x):= \frac{1}{2\pi} \int_{\partial \Omega} \ln |x-y| \; \mu (y) \, d\sigma(y), \quad x \in \mathbb{R}^2.
\end{equation}
The following formulas for the single layer potential are well known (see, for example, \cite{AmKa07Book2, Fo95}) :
\begin{align}
\mathcal{S}_{\partial \Omega} [\mu] \big|_+ (x) &= \mathcal{S}_{\partial \Omega} [\mu] \big|_- (x) \quad \text{ on } \partial \Omega, \label{singlejump1} \\
\partial_\nu \mathcal{S}_{\partial \Omega} [\mu] \big|_\pm (x) &= \biggl( \pm \frac{1}{2} I + \mathcal{K}_{\partial \Omega}^* \biggr) [\mu] (x) \quad \text{ on } \partial \Omega. \label{singlejump2}
\end{align}
Here the operator $\mathcal{K}_{\partial \Omega}^*$ is defined by
\begin{equation}\label{introkd}
\mathcal{K}_{\partial \Omega}^* [\mu] (x):=
\frac{1}{2\pi}\int_{\partial \Omega} \frac{(x-y) \cdot \nu_x}{|x-y|^2} \; \mu(y)\,d\sigma(y), \quad x \in \partial  \Omega,
\end{equation}
and is called the NP operator on $\partial \Omega$.

The NP operator $\mathcal{K}_{\partial \Omega}^*$ is not self-adjoint on $L^2(\partial\Omega)$ unless $\partial \Omega$ is a circle. However it can be realized as a self-adjoint operator on $H^{-1/2}(\partial\Omega)$ using the Plemelj's symmetrization principle:
\begin{equation}
\mathcal{S}_{\partial \Omega} \mathcal{K}_{\partial \Omega}^* = \mathcal{K}_{\partial \Omega} \mathcal{S}_{\partial \Omega}
\end{equation}
(see \cite{KPS}). The NP operator $\mathcal{K}^*_{\partial \Omega}$ is compact if $\partial \Omega $ is $C^{1,\alpha}$, and hence its eigenvalues are real, of finite multiplicities, and accumulate to 0.

We now consider the single layer potential on $\partial D= \partial D_1 \cup \partial D_2$. For this section we assume that $D_1$ and $D_2$ are bounded domains with Lipschitz continuous boundaries for some $\alpha >0$. For $\varphi := (\varphi_1, \varphi_2)^T \in \mathcal{H}^{-1/2}(\partial D)$, let
$$
\mathcal{S}_{\partial D}[\varphi](x) := \mathcal{S}_{\partial D_1}[\varphi_1](x) +  \mathcal{S}_{\partial D_2}[\varphi_2](x), \quad x \in \mathbb{R}^2.
$$
We seek a solution to \eqnref{mainproblem} (assuming that $\eta =(\eta_1,\eta_2)^T \in \mathcal{H}_0^{-1/2}(\partial D)$) in the form of
\begin{equation}\label{inteqn1}
u(x) := \mathcal{S}_{\partial D}[\varphi](x)
\end{equation}
for some $\varphi \in \mathcal{H}_0^{-1/2}(\partial D)$.

It follows from the jump relation \eqnref{singlejump2} and the jump condition of the flux (the third line in \eqnref{mainproblem}) that $\varphi$ should satisfy the following integral equations on $\partial D_1$ and $\partial D_2$, respectively:
$$
(1/2I+\mathcal{K}_{\partial D_1}^*)[\varphi_1] + \partial_{\nu} \mathcal{S}_{\partial D_2}[\varphi_2]
- k_1 \big( (-1/2I+\mathcal{K}_{\partial D_1}^*)[\varphi_1] + \partial_\nu \mathcal{S}_{\partial D_2}[\varphi_2] \big) = (k_1-1)\eta_1,
$$
and
$$
\partial_\nu \mathcal{S}_{\partial D_1}[\varphi_1] + (1/2I+ \mathcal{K}_{\partial D_2}^*)[\varphi_2] - k_2 \big(
\partial_{\nu} \mathcal{S}_{\partial D_1}[\varphi_1] + (-1/2I+ \mathcal{K}_{\partial D_2}^*)[\varphi_2] \big) = (k_2-1)\eta_2.
$$
These integral equations can be expressed in a short form as follows:
\begin{equation}\label{Gvfg}
\left( \Lambda - \mathbb{K}_{\partial D}^* \right) [\varphi] = \eta,
\end{equation}
where
\begin{equation}
\Lambda:= \begin{bmatrix} \lambda_1 I & 0 \\ 0 & \lambda_2 I \end{bmatrix}
\end{equation}
($\lambda_j$ is the constant defined by \eqnref{Glj} and $I$ is the identity operator), and
\begin{equation}
\mathbb{K}_{\partial D}^* \begin{bmatrix} \varphi_1 \\ \varphi_2 \end{bmatrix} = \begin{bmatrix} \mathcal{K}_{\partial D_1}^*[\varphi_1] & \partial_{\nu} \mathcal{S}_{\partial D_2}[\varphi_2]|_{\partial D_1} \\
\partial_{\nu} \mathcal{S}_{\partial D_1}[\varphi_1]|_{\partial D_2}  & \mathcal{K}_{\partial D_2}^*[\varphi_2] \end{bmatrix}.
\end{equation}
The operator $\mathbb{K}_{\partial D}^*$ is the NP operator associated with $\partial D=\partial D_1 \cup \partial D_2$.

Let, for $i,j=1,2$, $S_{ij}$ be the operator defined by
\begin{equation}
S_{ij} [\mu]:= \mathcal{S}_{\partial D_j}[\mu] |_{\partial D_i}, \quad \mu \in H^{-1/2}(\partial D_j),
\end{equation}
where the integration is carried out over $\partial D_j$ and the integral is evaluated on $\partial D_i$. Then $S_{ij}$ maps $H^{-1/2} (\partial D_j)$ into $H^{1/2} (\partial D_i)$. We then define
\begin{equation}
\mathbb{S}_{\partial D}: = \begin{bmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{bmatrix} .
\end{equation}
Then $\mathbb{S}_{\partial D}$ is a bounded operator from $\mathcal{H}^{-1/2}$ into $\mathcal{H}^{1/2}:= H^{1/2} (\partial D_1) \times H^{1/2} (\partial D_2)$.

The following analogy of the Plemelj's symmetrization principle was proved in \cite{ACKLM13}:
\begin{equation}\label{Ple}
\mathbb{S}_{\partial D} \mathbb{K}_{\partial D}^* = \mathbb{K}_{\partial D} \mathbb{S}_{\partial D}.
\end{equation}
We define a bilinear form on $\mathcal{H}^{-1/2}$ by
\begin{equation}\label{biform}
\langle \varphi, \psi \rangle_{\partial D}:= - (\varphi, \mathbb{S}_{\partial D}[\psi]),
\end{equation}
where the right-hand side is the duality pairing between $\mathcal{H}^{-1/2}$ and $\mathcal{H}^{1/2}$. Then $\langle \cdot, \cdot \rangle_{\partial D}$ is actually an inner product on $\mathcal{H}_0^{-1/2}$ and $\mathbb{K}_{\partial D}^*$ is compact and self-adjoint on $\mathcal{H}^{-1/2}$, namely,
$$
\langle \varphi, \mathbb{K}_{\partial D}^*[\psi] \rangle_{\partial D} = \langle \mathbb{K}_{\partial D}^*[\varphi], \psi \rangle_{\partial D},
$$
which is due to \eqnref{Ple}.

Then we have the following lemma. It is convenient to use the following notation: for a function $u$ on $\mathbb{R}^2$, $\partial_\nu u|_{\partial D_i \pm}$ denotes the trace of $\partial_\nu u$ on $\partial D_i$ from outside and inside of $D_i$, respectively, and
\begin{equation}
\partial u|_{\pm} := \begin{bmatrix} \partial_\nu u|_{\partial D_1 \pm} \\
\partial_\nu u|_{\partial D_2 \pm} \end{bmatrix}.
\end{equation}
If $u(x)= \mathcal{S}_{\partial D}[\varphi](x)$ for some $\varphi \in \mathcal{H}^{-1/2}$, then
\begin{equation}
\partial u|_\pm = (\pm 1/2I + \mathbb{K}_{\partial D}^*)[\varphi].
\end{equation}

\begin{lemma}\label{lem:1/2}
Let
\begin{equation}
E_{1/2}:= \{ \psi \in \mathcal{H}^{-1/2} : \langle \psi, \varphi \rangle_{\partial D}=0 \text{ for all $\varphi \in \mathcal{H}_0^{-1/2}$} \},
\end{equation}
so that $\mathcal{H}^{-1/2}= \mathcal{H}_0^{-1/2} \oplus E_{1/2}$.
\begin{itemize}
\item[(i)] Eigenvalues of $\mathbb{K}_{\partial D}^*$ on $\mathcal{H}_0^{-1/2}$ belongs to $(-1/2,1/2)$.
\item[(ii)] $\Lambda - \mathbb{K}_{\partial D}^*$ is invertible on $\mathcal{H}_0^{-1/2}(\partial D)$.
\item[(iii)] $E_{1/2}$ is the eigenspace of $\mathbb{K}_{\partial D}^*$ corresponding to the eigenvalue $1/2$ and consists of $\psi \in \mathcal{H}^{-1/2}$ such that $\mathcal{S}_{\partial D}[\psi]$ is constant on $D_1$ and $D_2$. In particular, $E_{1/2}$ is of two dimensions.
\end{itemize}
\end{lemma}
\begin{proof}
(i) is well-known (see, for example, \cite[Chapter XI, Section 11]{Kellog-book}).

Since $\mathbb{K}_{\partial D}^*$ is compact, (ii) follows from the injectivity of $\Lambda - \mathbb{K}_{\partial D}^*$ on $\mathcal{H}_0^{-1/2}(\partial D)$ by Fredholm alternative. If $\left( \Lambda - \mathbb{K}_{\partial D}^* \right) [\varphi]=0$, then $u(x) := \mathcal{S}_{\partial D}[\varphi](x)$ is the solution to \eqnref{mainproblem} with $\eta=0$ and $u(x)=O(|x|^{-1})$ as $|x| \to \infty$. It implies that $u \equiv 0$ and hence $\varphi=0$.

If $\psi \in E_{1/2}$, then $(\varphi, \mathbb{S}_{\partial D}[\psi])=0$ for all $\varphi \in \mathcal{H}_0^{-1/2}$. Thus $\mathbb{S}_{\partial D}[\psi]$ is constant, that is, $u:= \mathcal{S}_{\partial D}[\psi]$ is constant on $D_1$ and $D_2$. Therefore, we have
$$
\left( -1/2 I + \mathbb{K}_{\partial D}^* \right) [\psi] = \partial u|_- =0 ,
$$
and hence $\mathbb{K}_{\partial D}^* [\psi] = 1/2 \psi$. The converse can be proved by reversing the argument.
\end{proof}

\subsection{$1/2$ eigenvectors}\label{subsec:1/2}

Let
$q^1$ be the solution to
\begin{equation}\label{qonedef}
\begin{cases}
\displaystyle \Delta q^1  = 0 \quad & \mbox{in } \mathbb{R}^2 \setminus (\partial D_1 \cup \partial D_2),\\
\displaystyle q^1 = \mbox{constant} \quad &\mbox{on } \partial D_i, \ i=1,2,  \\
\noalign{\smallskip}
\displaystyle \int_{\partial D_i} \partial_\nu q^1|_+ \, d\sigma = (-1)^{i+1} , \quad &i=1,2, \\
\noalign{\smallskip}
\displaystyle q^1(x) =O(|x|^{-1}) \quad &\mbox{as } |x|\rightarrow \infty.
\end{cases}
\end{equation}
The existence of the solution $q^1$ (and $q^2$ below) is proved in \cite{ACKLY13} for inclusions $D_j$ of general shape, not necessarily disks.  Note that the constant values $q^1|_{\partial D_j}$ satisfy that $q^1|_{\partial D_1} <0$ and $q^1|_{\partial D_2} >0$, which can be shown using the maximum principle and Hopf's lemma. In particular, $q^1$ has a potential gap, namely,
\begin{equation}\label{qgap}
q^1|_{\partial D_2} - q^1|_{\partial D_1} >0.
\end{equation}
Let
\begin{equation}\label{mudef}
m= -\frac{q^1|_{\partial D_2}}{q^1|_{\partial D_1}} (>0),
\end{equation}
and let $q^2$ be the solution to
\begin{equation}\label{qtwodef}
\begin{cases}
\displaystyle \Delta q^2  = 0 \quad \mbox{in } \mathbb{R}^2 \setminus (\partial D_1 \cup \partial D_2),\\
\displaystyle q^2 = \mbox{constant} \quad \mbox{on } \partial D_i, \ i=1,2,  \\
\noalign{\smallskip}
\displaystyle \int_{\partial D_1} \partial_\nu q^2|_+ \, d\sigma = m, \quad \int_{\partial D_2} \partial_\nu q^2|_+ \, d\sigma = 1, \\
\displaystyle q^2 (x)= \frac{m+1}{2\pi} \ln|x| + O(|x|^{-1}) \quad\mbox{as }|x| \to \infty.
\end{cases}
\end{equation}

Let
\begin{equation}
\psi^i:=
\begin{bmatrix}
\partial_\nu q^i|_{\partial D_1+} \\ \partial_\nu q^i|_{\partial D_2+}
\end{bmatrix}, \quad j=1,2.
\end{equation}
One can see from \eqnref{qonedef} and \eqnref{qtwodef} that for $i=1,2$
$$
q^i(x)= \mathcal{S}_{\partial D}[\psi^i](x)+ \mathcal{D}_{\partial D}[q^i](x)
$$
for $x \in \mathbb{R}^2 \setminus (D_1 \cup D_2)$ by applying Green's theorem, where $\mathcal{D}_{\partial D}$ is the double layer potential, namely,
$$
\mathcal{D}_{\partial D} [\varphi] (x):=
\frac{1}{2\pi}\int_{\partial D} \frac{(x-y) \cdot \nu_x}{|x-y|^2} \; \varphi(y)\,d\sigma(y), \quad x \in \mathbb{R}^2 \setminus (\partial D_1 \cup \partial D_2).
$$
Since $q^i$ is constant on $\partial D_j$ for $j=1,2$, $\mathcal{D}_{\partial D}[q^i](x)=0$, and hence $q^i(x)= \mathcal{S}_{\partial D}[\psi^i](x)$ for $x \in \mathbb{R}^2 \setminus (D_1 \cup D_2)$. We then infer from Lemma \ref{lem:1/2} that $\psi^i$ belongs to the eigenspace $E_{1/2}$.

The eigenfunctions $\psi^1$ and $\psi^2$ constitute an orthogonal basis of $E_{1/2}$.
In fact, we have
\begin{align*}
\langle \psi^2 , \psi^1 \rangle_{\partial D} &= -(\psi^2, \mathbb{S}_{\partial D}[\psi^1]) \\
&= - q^1|_{\partial D_1} \int_{\partial D_1} \psi^2_1 d\sigma -q^1|_{\partial D_2} \int_{\partial D_2} \psi^2_2 d\sigma = -m q^1|_{\partial D_1} -  q^1|_{\partial D_2}=0,
\end{align*}
where the last equality follows from the definition \eqnref{mudef} of $m$.
Note that
$$
0= (\psi^1, \mathbb{S}_{\partial D}[\psi^2]) = q^2|_{\partial D_1} \int_{\partial D_1} \psi^1_1 \, d\sigma + q^2|_{\partial D_2} \int_{\partial D_2} \psi^1_2 \, d\sigma = q^2|_{\partial D_1} -  q^2|_{\partial D_2}.
$$
Thus we have
\begin{equation}\label{nogap}
q^2|_{\partial D_1} = q^2|_{\partial D_2},
\end{equation}
that is, $q^2$ has no potential gap on $\partial D_1$ and $\partial D_2$, from which one can infer that $\nabla q^2$ is bounded in $\mathbb{R}^2 \setminus (D_1 \cup D_2)$ (this can be proved by following the same lines of the proof of Theorem 2.1 in \cite{KLY}). It is also known that $|\nabla q^1| \gtrsim \epsilon^{-1/2}$ (see \cite{Yun} or the end of this subsection).

One can easily see that
\begin{equation}\label{psionenorm}
\| \psi^1 \|_{\partial D}^2 = q^1|_{\partial D_2} - q^1|_{\partial D_1},
\end{equation}
and
\begin{equation}\label{psitwonorm}
\| \psi^2 \|_{\partial D}^2 =  - q^2|_{\partial D_1} (1+ m)= q^2|_{\partial D_1} \left( \frac{q^1|_{\partial D_2}}{q^1|_{\partial D_1}} -1 \right).
\end{equation}
By \eqnref{qgap}, $\| \psi^1 \|_{\partial D} >0$. But we do not know if $q^2|_{\partial D_1} \neq 0$. If it is the case, then $\langle \cdot, \cdot \rangle_{\partial D}$ is an inner product not only on $\mathcal{H}_0^{-1/2}$ but also on $\mathcal{H}^{-1/2}$.

What we wrote so far in this subsection holds for inclusions of arbitrary shape (with Lipschitz continuous boundaries). If $D_1$ and $D_2$ are disks, then $q^1$ is explicitly given by
\begin{equation}\label{qBdef}
q^1(x) = \frac{1}{2\pi} (\ln|x-p_1|-\ln|x-p_2|), \quad x \in \mathbb{R}^2 \setminus (D_1 \cup D_2),
\end{equation}
where $p_1$ and $p_2$ be the unique fixed points of the combined reflections $R_1\circ R_2$ and $R_2\circ R_1$, respectively. Here, $R_j$ be the inversion with respect to the circle $\partial D_j$. Since $\partial D_1$ and $\partial D_2$ are the Apollonius circles with respect to points $p_1$ and $p_2$, $q^1$ is constant on $\partial D_1$ and $\partial D_2$ (constants may be different depending on the radii of the circles). This function was used for the first time in \cite{Yun} for analysis of field concentration. From the formula \eqnref{qBdef}, one can see that $\nabla q^1$ blows up at the order of $\epsilon^{-1/2}$ as $\epsilon$ tends to zero.

\subsection{Well-posedness for the transmission problems}

Let
\begin{equation}\label{Ximatrix}
\Xi := \begin{bmatrix} k_1-1 & 0 \\ 0 & k_2-1 \end{bmatrix}.
\end{equation}
For $\eta = (\eta_1, \eta_2)^T \in \mathcal{H}^{-1/2}(\partial D)$, let
$$
m_i := \frac{\langle \Xi \eta, \psi^i \rangle_{\partial D}}{\| \psi^i \|_{\partial D}^2} , \quad i=1,2,
$$
assuming that $\| \psi^2 \|_{\partial D} \neq 0$. One can see from \eqnref{psionenorm} that
\begin{equation}\label{Monedef}
m_1= \frac{- \sum_{j=1}^2 q^1|_{\partial D_j} (k_j-1) \int_{\partial D_j} \eta_j}{q^1|_{\partial D_2} - q^1|_{\partial D_1}},
\end{equation}
and from \eqnref{nogap} and \eqnref{psitwonorm} that
\begin{equation}\label{Mtwodef}
m_2= \frac{q^1|_{\partial D_1} \sum_{j=1}^2 (k_j-1) \int_{\partial D_j} \eta_j}{q^1|_{\partial D_2} - q^1|_{\partial D_1}}.
\end{equation}
We take \eqnref{Mtwodef} as the definition of $m_2$ when $\| \psi^2 \|_{\partial D} = 0$, or equivalently $q^2|_{\partial D_1} = 0$. We then have
\begin{equation}\label{E1/2compMj}
\Xi \eta - (m_1 \psi^1 + m_2 \psi^2) \in \mathcal{H}^{-1/2}_0(\partial D).
\end{equation}

Let $u$ be the solution to \eqnref{mainproblem} and let $v:= u- (m_1 q^1 + m_2 q^2)$. Then $v$ is the solution to \eqnref{mainproblem} with $\eta$ replace by $\eta^0$, which is defined by
\begin{equation}\label{gzerodef}
\eta^0:= \eta- \Xi^{-1} (m_1 \psi^1 + m_2 \psi^2).
\end{equation}
We emphasize that $\eta^0$ belongs to $\mathcal{H}_0^{-1/2}(\partial D)$.

\begin{prop}\label{prop:unique}
For $\eta \in \mathcal{H}^{-1/2}(\partial D)$, the solution $v$ to \eqnref{mainproblem} is unique and given by
\begin{equation}\label{mainprobsol}
v(x)= m_1 q^1(x) + m_2 q^2(x) + \mathcal{S}_{\partial D}[\varphi^0](x), \quad x \in \mathbb{R}^2,
\end{equation}
where $m_1, m_2$ are numbers given by \eqnref{Monedef} and \eqnref{Mtwodef}, and $\varphi^0 \in \mathcal{H}_0^{-1/2}(\partial D)$ is the unique solution to
\begin{equation}\label{Gvfgzero}
\left( \Lambda - \mathbb{K}_{\partial D}^* \right) [\varphi^0] = \eta^0
\end{equation}
with $\eta^0$ given by \eqnref{gzerodef}.
\end{prop}

\begin{proof}
It suffices to show that the solution is unique. For that, we prove that the solution $u$ to \eqnref{mainproblem} is zero assuming $\eta=0$. Let $v$ be the solution to \eqnref{mainproblem} with $\eta=0$. If $v (x)= c \ln |x| +O(|x|^{-1})$ for some constant $c$ as $|x| \to \infty$, let
$$
w(x):= v(x)- \frac{2\pi c}{m+1}q^2(x).
$$
Then, $w$ is a solution to \eqnref{mainproblem} with $\eta= -\frac{2\pi c}{m+1} \left(\frac{1}{k_1-1} \partial_\nu q^2|_{\partial D_1+}, \frac{1}{k_2-1} \partial_\nu q^2|_{\partial D_2+}\right)^T$ and $w(x)= O(|x|^{-1})$ as $|x| \to \infty$ by the fourth line of \eqnref{qtwodef}. It then follows from the third line of \eqnref{qtwodef} that
$$
\sum_{j=1}^2 \int_{\partial D_j} (\partial_\nu w|_+ - k_j \partial_\nu w|_-) = -\frac{2\pi c}{m+1} \sum_{j=1}^2 \int_{\partial D_j} \partial_\nu q^2|_{+} = -2\pi c.
$$
On the other hand, we have from divergence theorem that
$$
\sum_{j=1}^2 \int_{\partial D_j} (\partial_\nu w|_+ - k_j \partial_\nu w|_-) = -\int_{\mathbb{R}^2 \setminus D} \Delta w - \sum_{j=1}^2 k_j \int_{D_j} \Delta w =0.
$$
Thus, we have $c=0$, and hence $v (x)= O(|x|^{-1})$ as $|x| \to \infty$ from which we conclude that $v=0$.
\end{proof}

We have the following corollary:
\begin{cor}\label{cor:repin}
Let $F, \eta^0, m_1, m_2$ be quantities respectively given by \eqnref{Fdef}, \eqnref{gzerodef}, \eqnref{Monedef}, and \eqnref{Mtwodef}, and let $\varphi^0 \in \mathcal{H}_0^{-1/2}(\partial D)$ be the solution to \eqnref{Gvfgzero}.
Then, the unique solution $u$ to \eqnref{NCE} is given by
\begin{equation}\label{Nmainprobsol}
u= F + m_1 q^1 + m_2 q^2 + \mathcal{S}_{\partial D}[\varphi^0].
\end{equation}
\end{cor}

If $\eta_j= \partial_\nu H|_{\partial D_j}$ for $j=1,2$ where $H$ is the harmonic function appearing in \eqnref{HCE}, then
$$
\int_{\partial D_j} \eta_j = \int_{\partial D_j} \partial_\nu H = 0.
$$
Thus, in this case, we have $m_1=m_2=0$.

\begin{cor}\label{cor:repho}
Let $\varphi \in \mathcal{H}_0^{-1/2}(\partial D)$ be the solution to \eqnref{Gvfg} with $\eta_j= \partial_\nu H|_{\partial D_j}$.
Then, the unique solution $u$ to \eqnref{HCE} is given by
\begin{equation}\label{Hmainprobsol}
u= H + \mathcal{S}_{\partial D}[\varphi].
\end{equation}
\end{cor}

\section{NP spectrum on two circles}\label{sec:twocircles}

From now on we suppose that $D_1$ and $D_2$ are disks of radii $r_1$ and $r_2$, respectively. In this section we recall a complete spectrum of the NP operator on $\mathcal{H}^{-1/2}_0$ using the NP spectrum on concentric disks and M\"obius invariance of the NP spectrum (see, for example, \cite{Sch1}). NP spectrum on two disks is derived in \cite{BT} and we follow the method there. However, we need to review the method in some detail since we need to prove M\"obius invariance of the inner product (Lemma \ref{single_star} (ii)), and to derive the solution formula for the integral equation \eqnref{Gvfg} (Lemma \ref{lem:intsol}).

\subsection{NP spectrum on concentric circles}

Let $D_j^*$, $j=1,2$ be the disk of radius $R_j$ centered at $0$, where $R_1 < R_2$, and let $\mathbb{K}_{\partial D^*}$ be the NP operator on $\partial D^*$.
The complete spectrum of $\mathbb{K}_{\partial D^*}$ on $\mathcal{H}_0^{-1/2}(\partial D^*)$ is known (see \cite{ACKLM13}) as we describe in the following lemma.

\begin{lemma}\label{starspec}
Let $\rho:= R_1/R_2$. The eigenvalues of $\mathbb{K}_{\partial D^*}$ on $\mathcal{H}_0^{-1/2}(\partial D^*)$ are $\pm \frac{1}{2} \rho^{|n|}$ and the corresponding eigenfunctions are $f^{|n|,\pm}$ and $f^{-|n|,\pm}$ where
\begin{equation}\label{stareigenfunction}
f^{n,\pm} =  \left( \frac{1}{R_1} e^{in\theta} , \pm \frac{1}{R_2} e^{in\theta}  \right)^T , \quad n= \pm 1, \pm 2, \ldots.
\end{equation}
\end{lemma}

Actually Lemma \ref{starspec} is a consequence of the following formula whose proof can be found in \cite{ACKLM13}: if $\Omega$ is the disk of radius $R$ (centered at $0$), then for $n \neq 0$
\begin{equation}\label{slayerdisk}
\mathcal{S}_{\partial \Omega} \left[\frac{1}{R} e^{in\theta} \right](r, \theta) = \begin{cases}
\displaystyle -\frac{1}{2|n|} \left(\frac{r}{R}\right)^{|n|} e^{in\theta} \quad &\text{if } r \leq R, \\
\displaystyle -\frac{1}{2|n|} \left(\frac{R}{r}\right)^{|n|} e^{in\theta} \quad &\text{if } r>R.
\end{cases}
\end{equation}
Using \eqnref{slayerdisk}, one can see that the following formula holds:
if $n >0$, then
\begin{equation}\label{single_conc_circ2}
\mathcal{S}_{\partial D^*} [f^{n, \pm}] (\zeta) =
\begin{cases}
\displaystyle -\frac{1}{2n} (R_1^{-n} \pm R_2^{-n}) \zeta^n \quad & \text{if } |\zeta| \le R_1 , \\
\noalign{\smallskip}
\displaystyle -\frac{1}{2n} (R_1^n \bar{\zeta}^{-n} \pm R_2^{-n} \zeta^n ) \quad & \text{if } R_1 < |\zeta| \leq R_2, \\
\noalign{\smallskip}
\displaystyle -\frac{1}{2n} (R_1^{n} \pm R_2^{n}) \bar{\zeta}^{-n} \quad & \text{if } R_2 < |\zeta|,
\end{cases}
\end{equation}
and if $n <0$, then
\begin{equation}\label{single_conc_circ3}
\mathcal{S}_{\partial D^*} [f^{n, \pm}] (\zeta) = \overline{\mathcal{S}_{\partial D^*} [f^{-n, \pm}] (\zeta)},
\end{equation}
where $\zeta=r e^{i\theta}$. These formulas play a crucial role in what follows.

\subsection{M\"obius invariance and NP spectrum on two circles}

By translating and rotating $D$ if necessary, we may assume that their centers are located at $(c_1,0)$ and $(c_2,0)$, where
\begin{equation}\label{conedef}
c_1 = \frac{r_2^2 - r_1^2 - (r_1+r_2+\epsilon)^2}{2(r_1 + r_2 + \epsilon)} - \frac{\beta}{2} , \quad c_2 = c_1 + r_1 + r_2 + \epsilon
\end{equation}
with
\begin{equation}\label{Gbdef}
\beta = \sqrt{\epsilon} \frac{\sqrt{(2r_1 + \epsilon)(2r_2 +\epsilon)(2r_1 + 2r_2 + \epsilon)}}{r_1 + r_2 + \epsilon} ,
\end{equation}
where $\epsilon$ is the separation distance between $D_1$ and $D_2$.  These numbers are chosen so that $\partial D_1$ and $\partial D_2$ are mapped onto two concentric circles by an inversion in a circle and a translation:
if we let
\begin{equation}\label{zstardef}
z^*=Tz := \frac{\beta}{z} + 1,
\end{equation}
then $T(\partial D_j)$ ($j=1,2$) is the circle of the radius $R_j$ centered at $0$, where $R_j$ is given by
\begin{equation}\label{Ronetwodef}
R_1^2 = 1+ \frac{\beta}{c_1} , \quad
R_2^2 = 1+ \frac{\beta}{c_2} .
\end{equation}
Let
\begin{equation}
D_1^* := T(D_1)= \{ |\zeta| < R_1 \}, \quad D_2^* := T(D_2)=  \{ |\zeta| > R_2 \}.
\end{equation}

One can see from \eqnref{Gbdef} that
\begin{equation}\label{Gbasym}
\beta= \frac{4}{r_*} \sqrt{\epsilon} + O(\epsilon),
\end{equation}
where $r_*$ is defined by \eqnref{notation2}.
Since $c_1= -r_1 + O(\sqrt{\epsilon})$ and $c_2= r_2 + O(\sqrt{\epsilon})$, we see from \eqnref{Ronetwodef} that
\begin{equation}\label{Rest}
R_1 = 1 - \frac{\beta}{2r_1} + O(\epsilon), \quad
R_2 = 1 + \frac{\beta}{2r_2} + O(\epsilon),
\end{equation}
which together with \eqnref{Gbasym} yields
\begin{equation}\label{Grasymp}
\rho := \frac{R_1}{R_2} = 1 - r^* \sqrt{\epsilon} + O(\epsilon),
\end{equation}
as $\epsilon \to 0$.

Under the change of variables \eqnref{zstardef} the line and area elements are respectively transformed as
\begin{equation}\label{cvarea}
ds(z^*)= \frac{\beta}{|z|^2} ds(z)
\end{equation}
and
\begin{equation}\label{cvarea2}
dA(z^*)= \frac{\beta^2}{|z|^4} dA(z).
\end{equation}

Define the transformation $U$ from $\mathcal{H}^{-1/2}(\partial D)$ to $\mathcal{H}^{-1/2}(\partial D^*)$ by
\begin{align} \label{phistar}
(U\varphi)(z^*) =\varphi^*(z^*) := \frac{|z|^2}{\beta} \varphi(z), \quad z^* \in \partial D^*.
\end{align}
We obtain the following lemma (a similar lemma has been proved in \cite{JK} for the case with a single inclusion).

\begin{lemma}\label{single_star}
\begin{itemize}
\item[{\rm (i)}] The following transformation formula holds for any $\varphi \in \mathcal{H}^{-1/2}_0(\partial D)$:
\begin{equation}\label{singletrans}
\mathcal{S}_{\partial D^*} [ \varphi^* ](z^*)  = \mathcal{S}_{\partial D} [\varphi](z) - \mathcal{S}_{\partial D} [\varphi](0), \quad z \in \mathbb{C}.
\end{equation}
\item [{\rm (ii)}] $U$ is a unitary transform from $\mathcal{H}^{-1/2}_0(\partial D)$ to $\mathcal{H}^{-1/2}_0(\partial D^*)$, namely,
\begin{equation}\label{unitary}
\langle \varphi^*, \psi^* \rangle_{\partial D^*} = \langle \varphi, \psi \rangle_{\partial D}
\end{equation}
for any $\varphi, \psi \in \mathcal{H}^{-1/2}_0(\partial D)$.
\item [{\rm (iii)}] $\mathbb{K}^*_{\partial D}[\varphi] = \lambda \varphi$ ($|\lambda| < 1/2$) if and only if $\mathbb{K}^*_{\partial D^*}[\varphi^*] = -\lambda \varphi^*$. In fact, we have
\begin{equation}\label{KUUK}
\mathbb{K}^*_{\partial D^*} \circ U= - U \circ \mathbb{K}^*_{\partial D}.
\end{equation}
\end{itemize}
\end{lemma}
\begin{proof}
From the definition \eqref{zstardef} of $z^*$, we have
\begin{align}
|z^*-w^*| = \frac{\beta}{|z||w|}|z-w|,
\end{align}
and hence
\begin{align} \label{gammastar}
\Gamma(z^* - w^*) = \Gamma(z - w)  - \Gamma(w) - \Gamma(z) + \Gamma(\beta).
\end{align}
Here $\Gamma(z)= \frac{1}{2\pi} \ln|z|$. Hence for any $\varphi \in \mathcal{H}^{-1/2}_0(\partial D)$, we have from \eqref{cvarea} and \eqref{gammastar} that
\begin{align*}
\mathcal{S}_{\partial D^*}[\varphi^*](z^*) &= \int_{\partial D^*} \Gamma(z^* - w^*) \, \varphi^*(w^*) \, d\sigma(w^*) \\
&=\int_{\partial D} \left( \Gamma(z - w)  - \Gamma(w) - \Gamma(z) + \Gamma(\beta)\right) \, \varphi(w) \, d\sigma(w) \\
&= \mathcal{S}_{\partial D}[\varphi](z) - \mathcal{S}_{\partial D}[\varphi](0),
\end{align*}
which proves \eqref{singletrans}.

That $U$ is a unitary transform follows from \eqref{singletrans} immediately. In fact, we have
\begin{align*}
\langle \varphi^*, \psi^* \rangle_{\partial D^*} &= -\int_{\partial D^*} \varphi^*(z^*) \, \overline{\mathcal{S}_{\partial D^*}[\psi^*](z^*)} \, d\sigma(z^*) \\
&= -\int_{\partial D} \varphi(z) \; \overline{\left(\mathcal{S}_{\partial D}[\psi](z) - \mathcal{S}_{\partial D}[\psi](0)\right)} \, d\sigma(z) \\
&= \langle \varphi, \psi \rangle_{\partial D},
\end{align*}
where the last equality holds since $\varphi \in \mathcal{H}^{-1/2}_0(\partial D)$.

Let $n_z$ and $n_{z^*}$ be the complexified outward unit normal vectors on $\partial D_j$ and $\partial D^*_j$, respectively, namely,
$$
n_z = \frac{z-c_j}{r_j} \quad\mbox{if } z \in \partial D_j,
$$
and
$$
n_{z^*} = (-1)^{j+1} \frac{z^*}{R_j^*} \quad\mbox{if } z^* \in \partial D^*_j, \quad j=1,2.
$$
Then the following relation holds:
\begin{align}\label{nustar}
n_{z^*} = \frac{\bar{z}}{z} n_z.
\end{align}
In fact, since $|T(z)|^2=R_j$ if $z \in \partial D_j$, we see that ${\bar\partial_z} |T(z)|^2= {T(z)} \overline{\partial_z T(z)}$ is normal to $\partial D_j$. Thus,
$$
\frac{{\bar\partial_z} |T(z)|^2}{|{\bar\partial_z} |T(z)|^2|} = - \frac{T(z) z}{R_j \bar{z}} = c n_z, \quad \mbox{for } z \in \partial D_j,
$$
where $c$ is either $1$ or $-1$. One can easily see that $c=-1$, and hence \eqnref{nustar} holds.

Note that the outward normal derivatives $\nu_{z^*}$ and $\nu_z$ are given by
$$
\partial_{\nu_{z^*}}= 2\Re (n_{z^*} \partial_{z^*}), \quad \partial_{\nu_z}= 2\Re (n_{z} \partial_{z}).
$$
We thus have from \eqref{zstardef} and \eqnref{nustar} that
\begin{align} \label{nderistar}
\partial_{\nu_{z^*}} = 2\Re (n_{z^*} \partial_{z^*}) = -\frac{|z|^2}{\beta} 2\Re (n_{z} \partial_{z}) = -\frac{|z|^2}{\beta} \partial_{\nu_z} .
\end{align}
Hence we have from \eqref{gammastar} and \eqref{nderistar} the following identity
\begin{align*}
\mathbb{K}^*_{\partial D^*}[\varphi^*](z^*) &= \int_{\partial D^*} \partial_{\nu_{z^*}} \Gamma(z^* - w^*) \, \varphi^*(w^*) \, d\sigma(w^*) \\
&= -\frac{|z|^2}{\beta} \int_{\partial D} \partial_{\nu_z} \left( \Gamma(z - w) - \Gamma(z) \right) \, \varphi(w) \, d\sigma(w) \\
&= -\frac{|z|^2}{\beta} \mathbb{K}^*_{\partial D}[\varphi](z)
\end{align*}
which proves \eqref{KUUK}. This completes the proof.
\end{proof}

Let $f^{n, \pm}$ be the eigenfunctions of $\mathbb{K}_{\partial D^*}^*$ as given in \eqnref{stareigenfunction} and let $\varphi^{n, \pm}$, $n \in \mathbb{Z} \setminus \{0\}$ be functions such that
\begin{equation}\label{Gvfeigen}
U(\varphi^{n, \pm})= f^{n, \pm}.
\end{equation}
We have the following corollary from Lemma \ref{starspec}. This corollary is known (see \cite{BT}).

\begin{cor}
The eigenvalues of $\mathbb{K}_{\partial D}^*$ are $\pm\frac{1}{2} \rho^{|n|}$ and $1/2$, and
\begin{equation}\label{eigenconcir}
\mathbb{K}_{\partial D}^*[\varphi^{n, \pm}] = \mp \frac{1}{2} \rho^{|n|} \varphi^{n, \pm}.
\end{equation}
\end{cor}

Because of \eqnref{Grasymp}, more and more eigenvalues are approaching to $\pm 1/2$ as $\epsilon$ tends to $0$. This causes blow-up of the gradient (and higher derivatives) when conductivities of inclusions tend to $\infty$ or $0$. The phenomenon that eigenvalues are approaching to $\pm 1/2$ as $\epsilon$ tends to $0$ occurs when inclusions are of general shapes, as proved in \cite{BT2}.

\begin{lemma}
For each $n \neq 0$,
\begin{equation}\label{fnnorm}
\Vert \varphi^{n,\pm} \Vert^2_{\partial D} = \Vert f^{n,\pm} \Vert^2_{\partial D^*} =\frac{2\pi}{|n|} \left( 1 \pm \rho^{|n|} \right).
\end{equation}
\end{lemma}
\begin{proof}
The first identity in \eqnref{fnnorm} holds since $U$ is unitary. According to \eqnref{single_conc_circ2}, we have
\begin{align*}
\Vert f^{n,+} \Vert^2_{\partial D^*} &= - (f_{n,+}, \mathbb{S}_{\partial D^*}[f^{n,+}]) \\
&= \frac{1}{2|n|R_1} \int_{\partial D_1^*} \left(  1 + \rho^{|n|}  \right) d\theta + \frac{1}{2|n|R_2} \int_{\partial D_2^*} \left(  \rho^{|n|} + 1 \right) d\theta ,
\end{align*}
from which \eqnref{fnnorm} follows.
\end{proof}

We observe from the explicit forms \eqnref{stareigenfunction} of eigenfunctions $f^{n, \pm}$ and the formulas (\eqnref{single_conc_circ2} and \eqnref{single_conc_circ3}) of their single layer potentials on concentric circles that they are not only orthogonal to each other, but also satisfy the following property: if $n \neq m$, then
\begin{equation} \label{sepaorth}
\int_{\partial D_j^*} f_j^{n, \pm} \overline{\mathcal{S}_{\partial D^*}[f^{m, \pm}]} = \int_{\partial D_j^*} f_j^{n, \pm} \overline{\mathcal{S}_{\partial D^*}[f^{m, \mp}]} =0,
\end{equation}
and
\begin{equation} \label{sepa}
\int_{\partial D_j^*} f_j^{n, \pm} \overline{\mathcal{S}_{\partial D^*}[f^{n, +}]} = (-1)^{j+1} \int_{\partial D_j^*} f_j^{n, \pm} \overline{\mathcal{S}_{\partial D^*}[f^{n, -}]} \neq 0 ,
\end{equation}
for $j=1,2$. This property enables us to solve the integral equation \eqnref{Gvfg} by means of spectral decomposition.

\begin{lemma}\label{lem:intsol}
If
\begin{equation}
\eta= \sum_{n \neq 0} \big( C_{n, +} \varphi^{n, +} + C_{n, -} \varphi^{n, -} \big),
\end{equation}
then the solution $\varphi$ to \eqnref{Gvfg} is given by
\begin{equation}\label{Gvfsolexp}
\varphi = \sum_{n \neq 0} ( a_{n, +} \varphi^{n, +} + a_{n, -} \varphi^{n, -}),
\end{equation}
where
\begin{align}
a_{n,+} &= \frac{2}{4\lambda_1\lambda_2 - \rho^{2|n|}} \big( (\lambda_1 + \lambda_2 - \rho^{|n|}) C_{n, +} + (\lambda_1-\lambda_2) C_{n, -} \big), \label{solan+} \\
a_{n,-} &= \frac{2}{4\lambda_1\lambda_2 - \rho^{2|n|}} \big( -(\lambda_1-\lambda_2) C_{n, +} +  (\lambda_1 + \lambda_2 + \rho^{|n|}) C_{n, -} \big). \label{solan-}
\end{align}
\end{lemma}

\begin{proof}
One can see from \eqnref{eigenconcir}, \eqnref{sepaorth} and \eqref{sepa} that the following holds:
\begin{align*}
(\lambda_1 + 1/2 \rho^{|n|}) a_{n, +} + (\lambda_1 - 1/2 \rho^{|n|}) a_{n, -} &=  C_{n, +} + C_{n, -} , \\
(\lambda_2 + 1/2 \rho^{|n|}) a_{n, +} - (\lambda_2 - 1/2 \rho^{|n|}) a_{n, -} &=  C_{n, +} - C_{n, -}.
\end{align*}
Then, \eqnref{solan+} and \eqnref{solan-} immediately follow.
\end{proof}

If $\lambda_1=\lambda_2=\lambda$, then \eqnref{solan+} and \eqnref{solan-} are simplified to
$$
a_{n,\pm}= \frac{2C_{n, \pm}}{2\lambda \pm \rho^{|n|}},
$$
as expected.

\section{Spectral representations of the solution}\label{sec:specrep}

In this section we derive representation formulas for solutions to the problems \eqnref{NCE} and \eqnref{HCE} in terms of NP spectrum.

\subsection{Vanishing integral condition}

We first deal with the case when the integrals of $f$ over $D_1$ and $D_2$ vanish, namely,
\begin{equation}\label{intzero}
\int_{D_1} f = \int_{D_2} f = 0.
\end{equation}
In this case, $\partial_{\nu} F \in \mathcal{H}^{-1/2}_0(\partial D_j)$ for $j=1,2$ because of the condition \eqnref{intzero}. In fact, since $\Delta F= k_j^{-1} f$ on $D_j$ ($j=1,2$), we have
$$
\int_{\partial D_j} \partial_\nu F ds = \frac{1}{k_j} \int_{D_j} f=0.
$$
In particular, the constants $m_1$ and $m_2$ given in \eqnref{Monedef} and \eqnref{Mtwodef} become $0$. It thus follows from \eqnref{Nmainprobsol} that
the solution to \eqnref{NCE} admits the following representation:
\begin{equation}\label{uFv}
u(x)= F(x) + \mathcal{S}_{\partial D} [\varphi](x)
\end{equation}
where $\varphi \in \mathcal{H}^{-1/2}_0(\partial D)$ is the solution to the integral equation
\begin{equation}\label{GvfpF}
\left( \Lambda - \mathbb{K}_{\partial D}^* \right) [\varphi] = \partial F.
\end{equation}

Since $\partial_{\nu} F \in \mathcal{H}^{-1/2}_0(\partial D_j)$ for $j=1,2$, the following Neumann boundary value problem admits a unique solution which we denote by $H_j$:
\begin{equation}\label{bvp}
\begin{cases}
\Delta H_j = 0 \quad &\text{in } D_j, \\
\partial_{\nu} H_j = \partial_{\nu} F \quad &\text{on } \partial D_j.
\end{cases}
\end{equation}
Let $h_j$ be the analytic function in $D_j^*$ such that $h_1(0)=0$, $\lim_{|\zeta| \to \infty} h_2(\zeta)=0$, and
\begin{equation}\label{hjdef}
H_j(z)= \Re (h_j \circ T)(z)  + C_j, \quad z \in D_j
\end{equation}
for some constant $C_j$,  where $T$ is the transformation defined in \eqnref{zstardef} and $\Re$ denotes the real part.

We look into the integral equation \eqnref{GvfpF}. We have from the Green's identity, the jump formula \eqref{singlejump2} and the relations \eqnref{cvarea}, \eqnref{Gvfeigen}, \eqnref{eigenconcir} that
\begin{align*}
\langle \partial F, \varphi^{n,\pm} \rangle_{\partial D} &= - \sum_{j=1}^2 \int_{\partial D_j} \partial_{\nu} F \; \overline{\mathcal{S}_{\partial D}[\varphi^{n,\pm}]} \\
&= - \sum_{j=1}^2 \int_{\partial D_j} H_j \; \overline{\partial_{\nu} \mathcal{S}_{\partial D}[\varphi^{n,\pm}]|_-} \\
&= \frac{1}{2} \left( 1 \pm \rho^{|n|} \right) \sum_{j=1}^2 \int_{\partial D_j} H_j \; \overline{\varphi_j^{n,\pm}} \\
&= \frac{1}{4} \left( 1 \pm \rho^{|n|} \right) \sum_{j=1}^2 \int_{\partial D^*_j} (h_j+ \overline{h_j})  \; \overline{f_j^{n,\pm}}.
\end{align*}
It then follows from \eqnref{fnnorm} that
\begin{equation}\label{Cnpmdef}
\frac{\langle \partial F, \varphi^{n, \pm} \rangle_{\partial D}}{\Vert \varphi^{n, \pm} \Vert_{\partial D}^2} = \frac{|n|}{8\pi} C_{n, \pm}
:= \frac{|n|}{8\pi} \sum_{j=1}^2 \int_{\partial D^*_j} (h_j+ \overline{h_j})  \; \overline{f_j^{n,\pm}}.
\end{equation}
Thus we have
\begin{equation}
\partial F= \sum_{n \neq 0} \frac{|n|}{8\pi} \big( C_{n, +} \varphi^{n, +} + C_{n, -} \varphi^{n, -} \big).
\end{equation}

By Lemma \ref{lem:intsol}, the solution $\varphi$ is given by
\begin{equation}\label{Gvfexp}
\varphi= \sum_{n \neq 0} ( a_{n, +} \varphi^{n, +} + a_{n, -} \varphi^{n, -}),
\end{equation}
where
\begin{align}
a_{n,+} &= \frac{|n|}{4\pi( 4\lambda_1\lambda_2 - \rho^{2|n|})} \big( (\lambda_1 + \lambda_2 - \rho^{|n|}) C_{n, +} + (\lambda_1-\lambda_2) C_{n, -} \big), \label{an+} \\
a_{n,-} &= \frac{|n|}{4\pi( 4\lambda_1\lambda_2 - \rho^{2|n|})} \big( -(\lambda_1-\lambda_2) C_{n, +} +  (\lambda_1 + \lambda_2 + \rho^{|n|}) C_{n, -} \big). \label{an-}
\end{align}
The transformed function $\varphi^*=U(\varphi)$ is given by
\begin{align*}
\varphi^* &= \sum_{n \neq 0} \frac{|n|}{4\pi( 4\lambda_1\lambda_2 - \rho^{2|n|})} \Big[ C_{n, +} \Big( (\lambda_1 + \lambda_2 - \rho^{|n|})f^{n, +} - (\lambda_1-\lambda_2)f^{n, -} \Big) \\
& \qquad + C_{n, -} \Big( (\lambda_1-\lambda_2) f^{n, +} + (\lambda_1 + \lambda_2 + \rho^{|n|}) f^{n, -} \Big) \Big] .
\end{align*}
Thus, we have
\begin{align*}
\mathcal{S}_{\partial D^*} [\varphi^*] &= \sum_{n \neq 0} \frac{|n|}{4\pi( 4\lambda_1\lambda_2 - \rho^{2|n|})} \Big[ C_{n, +} \Big( (\lambda_1 + \lambda_2 - \rho^{|n|})H_{n, +} - (\lambda_1-\lambda_2) H_{n, -} \Big) \\
& \qquad + C_{n, -} \Big( (\lambda_1-\lambda_2) H_{n, +} + (\lambda_1 + \lambda_2 + \rho^{|n|}) H_{n, -} \Big) \Big] .
\end{align*}
Here we use notation $H_{n, \pm}(\zeta):= \mathcal{S}_{\partial D^*} [f^{n, \pm}](\zeta)$ for simplicity of expressions. Here and afterwards, we use $\zeta$ for $z^*=T(z)$.

Observe from the definition \eqnref{Cnpmdef} that $\overline{C_{n, \pm}}= C_{-n, \pm}$, and from \eqnref{single_conc_circ3} that $\overline{H_{n, \pm}}= H_{-n, \pm}$. We then have
\begin{equation}\label{SV1V2}
\mathcal{S}_{\partial D^*} [\varphi^*](\zeta) = \frac{V(\zeta) + \overline{V(\zeta)}}{2} ,
\end{equation}
where
\begin{align*}
V(\zeta) &= \sum_{n=1}^\infty \frac{n}{2\pi( 4\lambda_1\lambda_2 - \rho^{2n})} \Big[ C_{n, +} \Big( (\lambda_1 + \lambda_2 - \rho^{n})H_{n, +}(\zeta) - (\lambda_1-\lambda_2) H_{n, -}(\zeta) \Big) \\
& \qquad + C_{n, -} \Big( (\lambda_1-\lambda_2) H_{n, +}(\zeta) + (\lambda_1 + \lambda_2 + \rho^{n}) H_{n, -}(\zeta) \Big) \Big] .
\end{align*}

Since $h_1$ is analytic in $|\zeta| < R_1$ and $h_1(0)=0$, it can be expressed as
\begin{equation}\label{honeexp}
h_1(\zeta) = \sum_{m=1}^{\infty} a_{1,m} \zeta^m , \quad |\zeta| \leq R_1
\end{equation}
for some coefficients $a_{1,m}$. The function $h_2$ can be expressed as
\begin{equation}\label{htwoexp}
h_2(\zeta) = \sum_{m=1}^{\infty} a_{2,m} \zeta^{-m} , \quad |\zeta| \ge R_2
\end{equation}
for some coefficients $a_{2,m}$. Then, if $n >0$, then one can see from the formula \eqnref{stareigenfunction} for $f^{n,\pm}$ that
$$
C_{n, \pm} = \sum_{j=1}^2 \int_{\partial D^*_j} (h_j +\overline{h_j}) \; \overline{f_j^{n,\pm}} = 2\pi (a_{1,n} R_1^n \pm a_{2,n} R_2^{-n}).
$$
Thus $V(\zeta)$ can be written as
\begin{equation}\label{VB12}
V(\zeta) = - \frac{1}{2} \sum_{n=1}^\infty \frac{1}{4\lambda_1\lambda_2 - \rho^{2n}} \left( B_{1,n} + B_{2,n} \right) ,
\end{equation}
where
\begin{align}
B_{1,n}(\zeta) &= a_{1,n} R_1^n \Big( (2\lambda_1 - \rho^{n}) (-2n H_{n, +}(\zeta)) + (2 \lambda_2 + \rho^{n}) (-2n H_{n, -}(\zeta)) \Big), \label{B1n} \\
B_{2,n}(\zeta) &= \overline{a_{2,n}} R_2^{-n} \Big( (2\lambda_2 - \rho^{n})(-2n H_{n, +}(\zeta)) - (2 \lambda_1 + \rho^{n}) (-2n H_{n, -}(\zeta)) \Big) . \label{B2n}
\end{align}

We see from \eqnref{single_conc_circ2} that
$$
-2n H_{n, \pm}(\zeta) =
\begin{cases}
\displaystyle (R_1^{-1} \zeta)^n \pm (R_2^{-1} \zeta)^n \quad & \text{if } |\zeta| \le R_1 , \\
\noalign{\smallskip}
\displaystyle (R_1^{-1} \bar{\Gz})^{-n} \pm (R_2^{-1} \zeta)^n \quad & \text{if } R_1 < |\zeta| \leq R_2, \\
\noalign{\smallskip}
\displaystyle (R_1^{-1} \bar{\Gz})^{-n} \pm (R_2^{-1} \bar{\Gz})^{-n} \quad & \text{if } R_2 < |\zeta|.
\end{cases}
$$
Thus we have
\begin{align}
&B_{1,n}(\zeta) = 2a_{1,n} \times \nonumber \\
& \begin{cases}
\displaystyle (\lambda_1+\lambda_2) \zeta^n + (\lambda_1-\lambda_2) (\rho\zeta)^n - (\rho^2 \zeta)^n,  \quad  & |\zeta| \le R_1 , \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) (R_1^2 \bar{\Gz}^{-1})^n + (\lambda_1-\lambda_2) (\rho\zeta)^n - (\rho^2 \zeta)^n , \quad  & R_1 < |\zeta| \leq R_2, \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) (R_1^2 \bar{\Gz}^{-1})^n + (\lambda_1-\lambda_2) (R_1R_2 \bar{\Gz}^{-1})^n - (R_1^2 \bar{\Gz}^{-1})^n , \quad & R_2 < |\zeta|,
\end{cases} \label{Bltwo}
\end{align}
and
\begin{align}
&B_{2,n}(\zeta) = 2 \overline{a_{2,n}} \times \nonumber \\
&
\begin{cases}
\displaystyle (\lambda_1+\lambda_2) (R_2^2 \zeta^{-1})^{-n} - (\lambda_1-\lambda_2) (R_1R_2 \zeta^{-1})^{-n} - (R_2^2 \zeta^{-1})^{-n}, \quad  & |\zeta| \le R_1 , \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) (R_2^2 \zeta^{-1})^{-n} - (\lambda_1-\lambda_2) (\rho^{-1} \bar{\Gz})^{-n} - (\rho^{-2} \bar{\Gz})^{-n}, \quad & R_1 < |\zeta| \leq R_2, \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) \bar{\Gz}^{-n} - (\lambda_1-\lambda_2) (\rho^{-1} \bar{\Gz})^{-n} - (\rho^{-2} \bar{\Gz})^{-n},  \quad  & R_2 < |\zeta|.
\end{cases}\label{Blthree}
\end{align}

Let, for $j=1,2$,
\begin{equation}
A_j(\zeta):=-\frac{1}{2} \sum_{n=1}^\infty \frac{B_{j,n}(\zeta)}{4\lambda_1\lambda_2 - \rho^{2n}} .
\end{equation}
Then we have from \eqnref{VB12} that
\begin{equation}\label{A1A2}
V(\zeta) = A_1(\zeta)+ A_2(\zeta).
\end{equation}
It follows from \eqnref{Bltwo} that
\begin{equation}\label{A1def}
A_1(\zeta)
= \begin{cases}
\displaystyle (\lambda_1+\lambda_2) w_1(\zeta) & \\
\qquad + (\lambda_1-\lambda_2) w_1(\rho\zeta) - w_1(\rho^2 \zeta), \quad & |\zeta| \le R_1 , \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) w_1(R_1^2 \bar{\Gz}^{-1}) & \\
\qquad + (\lambda_1-\lambda_2) w_1(\rho\zeta) - w_1(\rho^2 \zeta), \quad & R_1 < |\zeta| \leq R_2, \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) w_1(R_1^2 \bar{\Gz}^{-1}) & \\
\qquad + (\lambda_1-\lambda_2) w_1(R_1R_2 \bar{\Gz}^{-1}) - w_1(R_1^2 \bar{\Gz}^{-1}),  \quad & R_2 < |\zeta|,
\end{cases}
\end{equation}
where
\begin{equation}\label{w1def}
w_1(\zeta) = \sum_{n=1}^\infty \frac{a_{1,n} \zeta^n}{4\lambda_1\lambda_2 - \rho^{2n}} , \quad |\zeta| \le R_1,
\end{equation}
and from \eqnref{Blthree} that
\begin{equation}\label{A2def}
A_{2}(\zeta) =
\begin{cases}
\displaystyle (\lambda_1+\lambda_2) w_2(R_2^2 \zeta^{-1}) & \\
\qquad - (\lambda_1-\lambda_2) w_2(R_1R_2 \zeta^{-1}) - w_2(R_2^2 \zeta^{-1}), \quad  & |\zeta| \le R_1 , \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) w_2(R_2^2 \zeta^{-1}) & \\
\qquad - (\lambda_1-\lambda_2) w_2(\rho^{-1} \bar{\Gz}) - w_2(\rho^{-2} \bar{\Gz}), \quad & R_1 < |\zeta| \leq R_2, \\
\noalign{\smallskip}
\displaystyle (\lambda_1+\lambda_2) w_2(\bar{\Gz}) & \\
\qquad - (\lambda_1-\lambda_2) w_2(\rho^{-1} \bar{\Gz}) - w_2(\rho^{-2} \bar{\Gz}),  \quad  & R_2 < |\zeta|,
\end{cases}
\end{equation}
where
\begin{equation}\label{w2def}
w_2(\zeta) = \sum_{n=1}^\infty \frac{a_{2,n} \zeta^{-n}}{4\lambda_1\lambda_2 - \rho^{2n}} , \quad |\zeta| \ge R_2.
\end{equation}

The functions $w_j$ can be expressed in terms of $h_j$. In fact, since
$$
\frac{1}{4\lambda_1\lambda_2 - \rho^{2n}} = \sum_{l=0}^\infty  \frac{\rho^{2ln}}{(4\lambda_1\lambda_2)^{l+1}},
$$
we have
\begin{equation}\label{w1def2}
w_1(\zeta) = \sum_{l=0}^\infty \frac{1}{(4\lambda_1\lambda_2)^{l+1}} \sum_{n=1}^\infty a_{1,n} (\rho^{2l}\zeta)^n = \sum_{l=0}^\infty \frac{h_1 (\rho^{2l}\zeta) }{(4\lambda_1\lambda_2)^{l+1}} , \quad |\zeta|<R_1.
\end{equation}
Likewise we have
\begin{equation}\label{w2def2}
w_2(\zeta) = \sum_{l=0}^\infty \frac{h_2 (\rho^{-2l}\zeta) }{(4\lambda_1\lambda_2)^{l+1}}, \quad |\zeta|> R_2.
\end{equation}

So far, we derived the representation formula for the solution as summarized in the following proposition.

\begin{prop}\label{prop:rep}
Suppose that $f$ satisfies \eqnref{intzero}. The solution $u$ to the inhomogeneous problem \eqnref{NCE} is represented as
\begin{equation}\label{firstrep}
u(z)= F(z) + \Re \left( A_{1}(T(z)) + A_2(T(z)) \right) + \mbox{const.},
\end{equation}
where $F, A_1, A_2$ are given by \eqnref{Fdef}, \eqnref{A1def}, \eqnref{A2def}, respectively, with $w_1$ and $w_2$ given by \eqnref{w1def2} and \eqnref{w2def2}.
\end{prop}

The solution $u$ to the homogeneous problem \eqnref{HCE} admits the same spectral representation formula as \eqnref{firstrep}. In this case, the solution $H_j$ of the boundary value problem \eqnref{bvp} is replaced with the given harmonic function $H$.

\subsection{General case}
We now deal with the more general case when $f$ does not necessarily satisfy the condition \eqnref{intzero}.

We still assume that the centers $c_1, \ c_2$ of $D_1, \ D_2$ lie on the real axis and satisfy $c_1<c_2$. Suppose that $c_1=0$ for convenience. Let
$$
U(r,\theta)=
\begin{cases}
\displaystyle \frac{k_1+1}{2k_1} r\cos \theta - \frac{(k_1-1)r_1^2}{2k_1} r^{-1} \cos \theta, \quad &r>r_1, \\
\noalign{\smallskip}
\displaystyle \frac{1}{k_1} r\cos\theta, \quad & r \le r_1,
\end{cases}
$$
which satisfies
\begin{equation}\label{Ueqn}
\begin{cases}
\Delta U = 0  \quad&\mbox{in } \mathbb{R}^2 \setminus \overline{D_1}, \\
\Delta U = 0  \quad&\mbox{in } D_1, \\
\displaystyle U|_+ - U|_- = 0  \quad&\mbox{on }\partial D_1 \\
\displaystyle \partial_\nu U|_+ - k_1 \partial_\nu U|_- = 0  \quad&\mbox{on }\partial D_1, \\
\displaystyle U (x) = \frac{k_1+1}{2k_1} r\cos \theta + O( |x|^{-1})~&\mbox{as }|x|\rightarrow\infty.
\end{cases}
\end{equation}
Let $\Theta(\theta)$ be a smooth function depending only on $\theta$. We assume that $\Theta(\theta)$ supported in $(\pi/2, 3\pi/2)$ and satisfies
\begin{equation}\label{3000}
r_1 \int_0^{2\pi} G(\theta) \cos \theta d\theta =1.
\end{equation}
Let $R(r)$ be a compactly supported smooth function of $r$ such that $R =1$ on $[0,r_1]$.
Let
$$
V_1(r, \theta):= \Theta(\theta)R(r) U(r,\theta),
$$
and
$$
v_1(r, \theta) := \begin{cases}
\nabla V_1, \quad &r>r_1, \\
k_1 \nabla V_1, \quad & r \le r_1.
\end{cases}
$$

\begin{lemma}\label{V1v1}
The function $V_1$ is supported in $\{ x_1 < c_1 \}$, is the solution to
\begin{equation}\label{V1eqn}
\begin{cases}
	\nabla \cdot \sigma \nabla V_1 = \nabla \cdot v_1 \quad &\text{in } \mathbb{R}^2, \\
	V_1(x) = O(|x|^{-1}) \quad &\text{as } |x| \to \infty,
\end{cases}
\end{equation}
and satisfies
\begin{equation}\label{V1est}
\| V_1 \|_{C^n(\overline{D_1})} \lesssim \frac{1}{k_1}, \quad \| V_1 \|_{C^n(\mathbb{R}^2 \setminus D_1)} \lesssim \frac{k_1+1}{k_1}
\end{equation}
for all positive integer $n$. Moreover, $v_1$ satisfies
\begin{equation}\label{v1est}
\| v_1 \|_{C^n(\overline{D_1})} \lesssim 1, \quad \| v_1 \|_{C^n(\mathbb{R}^2 \setminus D_1)} \lesssim \frac{k_1+1}{k_1}
\end{equation}
for all positive integer $n$ and
\begin{equation}\label{v1int}
\int_{D_1} \nabla \cdot v_1=1.
\end{equation}
\end{lemma}
\begin{proof}
It follows from the third and fourth lines of \eqnref{Ueqn} that
$$
V_1|_+ - V_1|_- = 0
$$
and
$$
\partial_\nu V_1|_+ - k_1 \partial_\nu V_1|_- = \Theta(\theta) (\partial_\nu U|_+ - k_1 \partial_\nu U|_- )= 0
$$
on $\partial D_1$. Since $V_1=0$ on $\{ x_1 \ge c_1 \}$, we have $V_1|_+ - V_1|_- = 0$ and $\partial_\nu V_1|_+ - k_2 \partial_\nu V_1|_- =0$ on $\partial D_2$. Thus $V_1$ satisfies \eqnref{V1eqn}. Estimates \eqnref{V1est} and \eqnref{v1est} immediately follow from the definitions of $V_1$ and $v_1$. For \eqnref{v1int}, we have
$$
\int_{D_1} \nabla \cdot v_1 = k_1 \int_{\partial D_1} \partial_\nu V_1 = \int_{\partial D_1} \Theta(\theta) \cos \theta d\sigma=1,
$$
where the last equality holds because of \eqnref{3000}.
\end{proof}

Likewise, we construct the functions $V_2$ such that $V_2$ and $v_2$, defined by
$$
v_2(r, \theta) := \begin{cases}
\nabla V_2, \quad &\mbox{in } D_2, \\
k_2 \nabla V_2, \quad &\mbox{in } \mathbb{R}^2 \setminus D_2,
\end{cases}
$$
satisfy the following lemma.

\begin{lemma}\label{V2v2}
The function $V_2$ is supported in $\{ x_1 > c_2 \}$, is the solution to
\begin{equation}\label{V2eqn}
\begin{cases}
	\nabla \cdot \sigma \nabla V_2 = \nabla \cdot v_2 \quad &\text{in } \mathbb{R}^2, \\
	V_2(x) = O(|x|^{-1}) \quad &\text{as } |x| \to \infty,
\end{cases}
\end{equation}
and satisfies
\begin{equation}\label{V2est}
\| V_2 \|_{C^n(\overline{D_2})} \lesssim \frac{1}{k_2}, \quad \| V_2 \|_{C^n(\mathbb{R}^2 \setminus D_2)} \lesssim \frac{k_2+1}{k_2}
\end{equation}
for all positive integer $n$. Moreover, $v_2$ satisfies
\begin{equation}\label{v2est}
\| v_2 \|_{C^n(\overline{D_2})} \lesssim 1, \quad \| v_2 \|_{C^n(\mathbb{R}^2 \setminus D_2)} \lesssim \frac{k_2+1}{k_2}
\end{equation}
for all positive integer $n$ and
\begin{equation}\label{v2int}
\int_{D_2} \nabla \cdot v_2=1.
\end{equation}
\end{lemma}

For a given $f=\nabla \cdot g$, let
\begin{equation}\label{fzero}
g_0:= g- \left( \int_{D_1} f \right) v_1 - \left( \int_{D_2} f \right) v_2 .
\end{equation}
Then $f_0=\nabla \cdot g_0$ satisfies \eqnref{intzero} and we arrive at the following proposition.

\begin{prop}\label{thm: decomp}
For a given $f=\nabla \cdot g$, let $f_0=\nabla \cdot g_0$ be defined by \eqnref{fzero}. The solution $u$ to \eqnref{NCE} can be decomposed as
\begin{equation}\label{usol2}
u= \left( \int_{D_1} f \right) V_1 + \left( \int_{D_2} f \right) V_2 + u_0,
\end{equation}
where $u_0$ is the solution to \eqnref{NCE} with $f=f_0$ which is of the form \eqnref{firstrep}.
\end{prop}

\section{Proofs of theorems}\label{sec:pfs}

In this subsection we prove Theorem \ref{thm:main1}-\ref{thm:main5}. Thanks to Proposition \ref{thm: decomp}, for the case of the problem \eqnref{NCE} we may assume that \eqnref{intzero} holds. In fact, since
$$
\left| \int_{D_j} f \right| = \left| \int_{\partial D_j} g \cdot \nu \right| \lesssim \| g \|_{L^\infty(\partial D_j)}
$$
for $j=1,2$, we have
$$
\| g_0 \|_{n,\alpha}^* \lesssim \| g \|_{n,\alpha}^*.
$$
Moreover, we infer from \eqnref{V1est}, \eqnref{V2est}, and \eqnref{usol2} that
$$
\| u \|_{n,0} \lesssim \frac{1}{k_1} \| g \|_{L^\infty(\partial D_1)} + \frac{1}{k_2} \| g \|_{L^\infty(\partial D_2)} + \| u_0 \|_{n,0}.
$$
Thus it is enough to prove Theorem \ref{thm:main1} and \ref{thm:main3} with $u$ and $g$ replaced with $u_0$ and $g_0$.

In the rest of this section, we assume that \eqnref{intzero} holds for the problem \eqnref{NCE}.
We begin by proving a preliminary lemma, for which we introduce the following transformations: With constants $\beta$ and $\rho$ appearing in \eqnref{conedef} and \eqnref{Grasymp}, let, for $\rho^2 \le t \le 1$ and $l=0,1,2, \ldots$,
\begin{equation}\label{Phidef}
\Phi_{l,t}(z):= T^{-1} (t \rho^{2l} T(z))= \frac{\beta z}{\beta t \rho^{2l} - (1-t\rho^{2l}) z}, \quad z \in D_1
\end{equation}
and
\begin{equation}\label{Psidef}
\Psi_{l,t}(z):= T^{-1}(t \rho^{-2l-1} T(z))= \frac{\beta z}{\beta t \rho^{-2l-1} - (1-t\rho^{-2l-1}) z}, \quad z \in D_2.
\end{equation}
Note that $\Phi_{l,t}$ maps $D_1$ into $D_1$ and $\Psi_{l,t}$ maps $D_2$ into $D_2$.

\begin{lemma}\label{lem:est1}
For $j=1,2$, let $V_j$ be an analytic function defined in $D_j$ such that $V_j \in C^{n}(\overline{D_j})$ for some positive integer $n$. For $\rho^2 \le t \le 1$ and $l=0,1,2, \ldots$, let $v_1(z) = (V_1 \circ \Phi_{l,t})(z)$ for $z \in D_1$ and $v_{2}(z) = (V_2 \circ \Psi_{l,t})(z)$ for $z \in D_2$. There is a constant $C$ independent of $t$ and $l$ such that
\begin{equation}\label{vjlest}
\| v_{j} \|_{C^{n}(\overline{D_j})} \le C \rho^{2l} (l+1)^{n} \| V_j \|_{C^{n}(\overline{D_j})}
\end{equation}
for $j=1,2$.
\end{lemma}

\begin{proof}
One can easily see from \eqnref{Phidef} that
\begin{equation}\label{Phideri}
\Phi_{l,t}^{(m)}(z) = \frac{m! \beta^2 t\rho^{2l} (1- t\rho^{2l})^{m-1}}{(\beta t\rho^{2l} - (1-t\rho^{2l}) z)^{m+1}}
\end{equation}
for $m=1,2, \ldots$, which can be written in terms of the $\zeta=T(z)$ variable as
\begin{equation}\label{Phideri2}
\Phi_{l,t}^{(m)}(z) = \beta^{1-m} t\rho^{2l} (1-t\rho^{2l})^{m-1} \frac{(\zeta-1)^{m+1}}{(t\rho^{2l} \zeta - 1)^{m+1}}.
\end{equation}
Since $|\zeta| \le R_1 <1$ and $0< t\rho^{2l}<1$, we have
\begin{equation}\label{4209}
|\zeta-1| \le 2 |t\rho^{2l} \zeta - 1 |,
\end{equation}
and hence
\begin{equation}\label{4210}
|\Phi_{l,t}^{(m)}(z)| \le 2^{m+1} \beta^{1-m} t\rho^{2l} (1-t\rho^{2l})^{m-1}.
\end{equation}
Since $\rho \le t \le 1$, we infer from \eqnref{Grasymp} that if $l \ge 1$
\begin{equation}\label{tGr2l}
1-t\rho^{2l} \sim 1-\rho^{2l+1} \sim (l+1) \sqrt{\epsilon}.
\end{equation}
If $l=0$, we have
\begin{equation}\label{tGr2l2}
1-t\rho^{2l} \lesssim \sqrt{\epsilon}.
\end{equation}
Since $\beta \approx \sqrt{\epsilon}$ (see \eqnref{Gbasym}), it follows from \eqnref{4210} that
\begin{equation}\label{Gest}
|\Phi_{l,t}^{(m)}(z)| \lesssim \rho^{2l} (l+1)^{m}
\end{equation}
for $m=0,1,2, \ldots$.

The desired estimate \eqnref{vjlest} for $j=1$ follows from \eqnref{Gest}. In fact,
we have
$$
v_{1}^{(n)}(z)= V_1^{(n)}(\Phi_{l,t}(z)) \Phi_{l,t}'(z)^n + \sum_{k=1}^{n-1} V_1^{(k)}(\Phi_{l,t}(z)) \sum_{\sum_{s=1}^t k_s m_s=n-k} C_{k,s} \Phi_{l,t}^{(k_s)}(z)^{m_s}
$$
for some constant $C_{k,s}$. It then follows from \eqnref{Gest} that
\begin{align*}
\| v_{1}^{(n)} \|_{C^{0}(\overline{D_1})} \lesssim \| V_1 \|_{C^{n}(\overline{D_1})} \left( \rho^{2ln} (l+1)^{n} + \sum_{k=1}^{n-1} \sum_{\sum_{s=1}^t k_s m_s=n-k} \rho^{2lm_s} (l+1)^{k_s m_s} \right),
\end{align*}
which leads us to \eqnref{vjlest} for $j=1$.

To prove \eqnref{vjlest} for $j=2$, we see as before that
\begin{equation}
\Psi_{l,t}^{(m)}(z) = \frac{m! \beta^2 t\rho^{-2l-1} (1- t\rho^{-2l-1})^{m-1}}{(\beta t\rho^{-2l-1} - (1-t\rho^{-2l-1}) z)^{m+1}}
\end{equation}
for $m=1,2, \ldots$, which can be written in terms of the $\zeta=T(z)$ variable as
$$
\Psi_{l,t}^{(m)}(z) = \beta^{1-m} t\rho^{-2l-1} (1-t\rho^{-2l-1})^{m-1} \frac{(\zeta-1)^{m+1}}{(t\rho^{-2l-1} \zeta - 1)^{m+1}}.
$$
Since $|\zeta| \ge R_2 >1$ and $1 \ge t\rho^{-2l-1}$, we have
$$
|\zeta-1| \le 2 |\zeta - t^{-1}\rho^{2l+1} |,
$$
and hence
$$
\left| \frac{(\zeta-1)^{m+1}}{(t\rho^{-2l-1} \zeta - 1)^{m+1}} \right| \le \rho^{2l} 2^{m+1}.
$$
As before, we have
\begin{equation}
t\rho^{-2l-1} -1  \sim (l+1)\beta.
\end{equation}
Thus we have
\begin{equation}\label{Psiest}
|\Psi_{l,t}^{(m)}(z)| \lesssim \rho^{2l} (l+1)^{m}
\end{equation}
for $m=0,1,2, \ldots$.

The estimate \eqnref{vjlest} for $j=2$ can be proved similarly to the case when $j=1$ using \eqnref{Psiest}.
\end{proof}

\medskip

\noindent{\sl Proof of Theorem \ref{thm:main1}}.
We first note that $4\lambda_1\lambda_2 >0$. Let $\tilde{h}_j := h_j \circ T$ so that $\Re \tilde{h}_j=H_j$ (see \eqnref{hjdef}). We see from \eqnref{w1def2} that
\begin{equation}\label{w1circT}
(w_1 \circ T)(z) = \sum_{l=0}^\infty \frac{\tilde{h}_1 (T^{-1} \rho^{2l} T(z)) }{(4\lambda_1\lambda_2)^{l+1}}.
\end{equation}
It follows from \eqnref{vjlest} that
$$
\| w_1 \circ T \|_{C^{n+1}(\overline{D_1})} \lesssim \sum_{l=0}^\infty
\frac{\rho^{2l} (l+1)^{n+1}}{(4\lambda_1\lambda_2)^{l+1}} \| H_1 \|_{C^{n+1}(\overline{D_1})} .
$$
Since
$$
\sum_{l=0}^\infty \frac{\rho^{2l} (l+1)^{n+1}}{(4\lambda_1\lambda_2)^{l+1}} \lesssim \frac{1}{(4\lambda_1\lambda_2 - \rho)^{n+1}},
$$
we have
\begin{equation}\label{w1Tz}
\| \nabla (w_1 \circ T) \|_{C^{n, \alpha}(\overline{D_1})} \lesssim \frac{\| H_1 \|_{C^{n+1}(\overline{D_1})}}{(4\lambda_1\lambda_2 - \rho)^{n+1}}.
\end{equation}
Likewise, we have
\begin{equation}\label{w2Tz}
\| \nabla (w_2 \circ T) \|_{C^{n, \alpha}(\overline{D_2})} \lesssim \frac{\| H_2 \|_{C^{n+1}(\overline{D_2})}}{(4\lambda_1\lambda_2 - \rho)^{n+1}}.
\end{equation}

One can show that all the terms in the expressions \eqnref{A1def} and \eqnref{A2def} of the functions $A_1$ and $A_2$ can be estimated using \eqnref{w1Tz} and \eqnref{w2Tz}. For example, the term $w(R_1^2 \bar{\Gz}^{-1})$ for $R_1 < |\zeta|$ (in \eqnref{A1def}) can be estimated as follows: Note that
$$
w_1(R_1^2 \bar{\Gz}^{-1}) = w_1(R_1^2 \overline{T(z)}^{-1}) = (w_1 \circ T) (G(z)),
$$
where
\begin{equation}\label{Vone}
G(z) := T^{-1} (R_1^2 \overline{T(z)^{-1}}).
\end{equation}
Note that
$$
G(z)= \frac{\beta(\beta+z)}{(R_1^2-1) z - \beta}.
$$
Thus, we have
$$
G^{(n)}(z)= -\frac{n! \beta^2 R_1^2 (R_1^2-1)^{n-1}}{(\beta -(R_1^2-1) z)^{n+1}}= - n!\beta^{1-n} R_1^2 (R_1^2-1)^{n-1}  \frac{(\zeta-1)^{n+1}}{(\zeta - R_1^2)^{n+1}}.
$$
Since $R_1 \le |\zeta|$ and $R_1<1$, we have
$$
|\zeta-1| \lesssim |\zeta - R_1^2|.
$$
Since $R_1=1+O(\sqrt{\epsilon})$ as shown in \eqnref{Rest}, we have
$$
|R_1^2-1| \lesssim \sqrt{\epsilon}.
$$
Thus we have
\begin{equation}\label{Voneest}
|G^{(n)}(z)| \lesssim 1
\end{equation}
for all $z \in \mathbb{R}^2 \setminus D_1$. Since $G$ maps $\mathbb{R}^2 \setminus D_1$ into $D_1$, we have from \eqnref{w1Tz} and \eqnref{Voneest} that
\begin{equation}
\| (w_1 \circ T) \circ G \|_{C^{n+1}(\mathbb{R}^2 \setminus D_1)} \lesssim \| w_1 \circ T \|_{C^{n+1}(\overline{D_1})} \lesssim \frac{\| H_1 \|_{C^{n+1}(\overline{D_1})}}{(4\lambda_1\lambda_2 - \rho)^{n+1}}.
\end{equation}

So far we showed that the following estimate holds independently of $k_j$ and $\epsilon$ for $j=1,2$:
$$
\| \nabla A_j \|_{n,0} \lesssim \frac{\| H_j \|_{C^{n+1}(\overline{D_1})}}{(4\lambda_1\lambda_2 - \rho)^{n+1}}.
$$
It follows from \eqnref{firstrep} that
\begin{equation}\label{6000}
\| \nabla (u-F) \|_{n,0} \lesssim \frac{\| H_1 \|_{C^{n+1}(\overline{D_1})} + \| H_2 \|_{C^{n+1}(\overline{D_2})}}{(4\lambda_1\lambda_2 - \rho)^{n+1}} .
\end{equation}

If $f=\nabla \cdot g$ and $g \in C^{n, \alpha}(\mathbb{R}^2)$ for some nonnegative integer $n$ and $0<\alpha<1$, then $F \in C^{n+1, \alpha}(\mathbb{R}^2)$ and
\begin{equation}\label{CZest}
\| F\|_{n+1, \alpha} \lesssim \| g \|_{n,\alpha}^*
\end{equation}
by the Calder\'on-Zygmund estimates. By regularity estimates of the boundary value problem, we have $H_j \in C^{n+1, \alpha}(\overline{D_j})$ and
\begin{equation}\label{H1est}
\| H_j \|_{C^{n+1, \alpha}(\overline{D_j})} \lesssim \| F\|_{C^{n+1, \alpha}(\overline{D_j})} \lesssim \| g \|_{n,\alpha}^*.
\end{equation}
It thus follows from \eqnref{6000} that
\begin{equation}\label{nablauest}
\| u \|_{n+1,0} \lesssim \frac{\| g \|_{n,\alpha}^*}{(4\lambda_1\lambda_2 - \rho)^{n+1}} .
\end{equation}
Since
\begin{equation}\label{Grest}
\rho \sim 1-r_* \sqrt{\epsilon}
\end{equation}
by \eqnref{Grasymp}, the desired estimate \eqnref{mainest1} follows from \eqnref{nablauest}.

A proof of optimality of \eqnref{mainest1} will be given in Section \ref{sec:opti}. This completes the proof of Theorem \ref{thm:main1}.  $\Box$

\bigskip

\noindent{\sl Proof of Theorem \ref{thm:main3}}.
If $(k_1-1)(k_2-1) <0$, then $\lambda_1\lambda_2 <0$. Let $\lambda:= -4\lambda_1\lambda_2$ and $Y_l^t(z):= \tilde{h}_1 (T^{-1} (t\rho^{2l} T(z))) = (\tilde{h}_1 \circ \Phi_{l,t})(z) $ ($\tilde{h}_j := h_j \circ T$ as before). We can rewrite \eqnref{w1circT} as
\begin{align*}
(w_1 \circ T)(z) &= -\frac{1}{\lambda} \sum_{l=0}^\infty \frac{(-1)^{l} Y_l^1(z)}{\lambda^{l}} \nonumber \\
&= -\frac{1}{\lambda}\sum_{i=0}^\infty \left( \frac{1}{\lambda^{2i}} - \frac{1}{\lambda^{2i+1}} \right) Y_{2i}^1(z) +  \frac{1}{\lambda^{2i+1}} (Y_{2i}^1(z)-Y_{2i+1}^1(z)) \\
&=: -\frac{1}{\lambda}\sum_{i=0}^\infty (I_i(z) + J_i(z)). \label{w1circT2}
\end{align*}

Since $\lambda >1$, we have
$$
\left| \frac{1}{\lambda^{2i}} - \frac{1}{\lambda^{2i+1}} \right| = \frac{\lambda-1}{\lambda^{2i+1}} \le \frac{\lambda-1}{\lambda^{i}}.
$$
It thus follows from \eqnref{vjlest} that
$$
\| I_i \|_{C^{n+1}(\overline{D_1})} \lesssim (\lambda-1) \frac{\rho^{i} (i+1)^{n+1}}{\lambda^{i}} \| H_1 \|_{C^{n+1}(\overline{D_1})}.
$$
Thus we have
\begin{equation}
\left\| \sum_{i=0}^\infty I_i \right\|_{C^{n+1}(\overline{D_1})} \lesssim
\frac{\lambda-1}{(\lambda-\rho)^{n+1}} \| H_1 \|_{C^{n+1}(\overline{D_1})},
\end{equation}
which together with \eqnref{H1est} yields
\begin{equation}\label{Isum}
\left\| \sum_{i=0}^\infty I_i \right\|_{C^{n+1}(\overline{D_1})} \lesssim
\frac{\| g \|_{n,\alpha}^*}{(\lambda-\rho)^{n}} .
\end{equation}

On the other hand, we have
\begin{align}
Y_{2i}^1(z)-Y_{2i+1}^1(z) &= Y_{2i}^1(z)-Y_{2i}^{\rho^2}(z) \nonumber \\
&= \int_{\rho^2}^1 \frac{\partial}{\partial t} Y_{2i}^t(z)dt =
\int_{\rho^2}^1 h_1' (\Phi_{2i,t}(z)) \frac{\partial}{\partial t} \Phi_{2i,t}(z) dt . \label{4220}
\end{align}
We see from \eqnref{Phideri2} that
\begin{align*}
\frac{\partial}{\partial t} \Phi_{2i,t}^{(m)}(z) &= \frac{\beta^{1-m} \rho^{4i} (\zeta-1)^{m+1}}{(t\rho^{4i} \zeta - 1)^{m+1}}  \\
& \quad \times  \left[ (1-t\rho^{4i})^{m-1} - (m-1) t\rho^{4i}(1-t\rho^{4i})^{m-2} - \frac{(m+1)\rho^{4i} \zeta}{t\rho^{4i} \zeta - 1} \right] \\
& =: M_1 + M_2 + M_3.
\end{align*}
Using \eqnref{4209} and \eqnref{tGr2l}, we see that
$$
|M_1| \lesssim \rho^{i} (i+1)^{m}, \quad
|M_2| \lesssim \frac{\rho^{i} (i+1)^{m-1}}{\sqrt{\epsilon}}.
$$
Since
$$
|t\rho^{4i} \zeta -1| \ge 1- t\rho^{i} R_1
$$
for any $\zeta$ with $|\zeta| \le R_1$, we have
$$
|t\rho^{4i} \zeta -1| \gtrsim \sqrt{\epsilon} (i+1).
$$
Thus we have
$$
|M_3| \lesssim \frac{\rho^{i} (i+1)^{m-1}}{\sqrt{\epsilon}}.
$$
It then follows that
\begin{equation}\label{4221}
\sum_{z \in D_1} \left| \frac{\partial}{\partial t} \Phi_{2i,t}^{(m)}(z) \right| \lesssim \rho^{i} (i+1)^{m} + \frac{\rho^{i} (i+1)^{m-1}}{\sqrt{\epsilon}}.
\end{equation}

By \eqnref{vjlest}, we have
$$
\| g \circ \Phi_{2i,t} \|_{C^{n+1}(\overline{D_1})} \lesssim \rho^{i} i^{n+1} \| g \|_{C^{n+1}(\overline{D_1})}
\lesssim \rho^{i} i^{n+1} \| g \|_{n+1,\alpha}^*.
$$
We then infer from \eqnref{4220} and \eqnref{4221} that
$$
\| Y_{2i}^1-Y_{2i+1}^1 \|_{C^{n+1}(\overline{D_1})} \lesssim (1-\rho^2) \left( \rho^{i} (i+1)^{m} + \frac{\rho^{i} (i+1)^{m-1}}{\sqrt{\epsilon}} \right) \| g \|_{n+1,\alpha}^*.
$$
Since $1-\rho^2 \lesssim \sqrt{\epsilon}$, we have
\begin{equation}
\| Y_{2i}^1-Y_{2i+1}^1 \|_{C^{n+1}(\overline{D_1})} \lesssim \rho^{i} \left( \sqrt{\epsilon} (i+1)^{n+1} + (i+1)^{n} \right) \| g \|_{n+1,\alpha}^*.
\end{equation}
It then follows that
$$
\left\| \sum_{i=0}^\infty J_i \right\|_{C^{n+1}(\overline{D_1})} \lesssim \sum_{i=0}^\infty \frac{\rho^{i}}{\lambda^i} \left( \sqrt{\epsilon} (i+1)^{n+1} + (i+1)^{n} \right) \| g \|_{n+1,\alpha}^* .
$$
Since
\begin{align*}
\sum_{i=0}^\infty \frac{\rho^{i}}{\lambda^i} \left( \sqrt{\epsilon} (i+1)^{n+1} + (i+1)^{n} \right)
&\lesssim \frac{\sqrt{\epsilon}}{(\lambda-\rho)^{n+1}} + \frac{1}{(\lambda-\rho)^{n}} \lesssim \frac{1}{(\lambda-\rho)^{n}},
\end{align*}
we have
\begin{equation}
\left\| \sum_{i=0}^\infty J_i \right\|_{C^{n+1}(\overline{D_1})} \lesssim \frac{\| g \|_{n+1,\alpha}^*}{(\lambda-\rho)^{n}}.
\end{equation}
This together with \eqnref{Isum} yields
\begin{equation}\label{4230}
\left\| w_1 \circ T \right\|_{C^{n+1}(\overline{D_1})} \lesssim \frac{\| g \|_{n+1,\alpha}^*}{(\lambda-\rho)^{n}}.
\end{equation}

Similarly, one can show that
\begin{equation}\label{4231}
\left\| w_2 \circ T \right\|_{C^{n+1}(\overline{D_2})} \lesssim \frac{\| g \|_{n+1,\alpha}^*}{(\lambda-\rho)^{n}}.
\end{equation}
Then, as before, one can show
$$
\| A_j \|_{n+1,0} \lesssim \frac{\| g \|_{n+1,\alpha}^*}{(\lambda-\rho)^{n}}
$$
for $j=1,2$, which yields the desired estimate \eqnref{mainest3-1}.
This completes the proof.  $\Box$

\bigskip

Since the solution to \eqnref{HCE} admits the same representation as the one in Proposition \ref{prop:rep}, Theorem \ref{thm:main4} and \ref{thm:main5} can be proved in the same way as Theorem \ref{thm:main1} and \ref{thm:main3}.

\section{Optimality}\label{sec:opti}

In this section we show that there is $f$ (or $H$) such that the reverse inequality holds in \eqnref{mainest1}, \eqnref{mainest4-1}, \eqnref{mainest3-1}, and \eqnref{mainest4-2}.

Let $F$ be a smooth function in $\mathbb{R}^2$ with a compact support such that $F(z)=x_1$ in a neighborhood of $\overline{D_1 \cup D_2}$. Let $f:= \Delta F$. Then the relation \eqnref{Fdef} holds, namely, $F$ is the Newtonian potential of $f$. Note that $f=0$ in $D_1 \cup D_2$. The solution $H_j$ of \eqnref{bvp} is given by $H_j(z)=x_1$ in $D_j$ for $j=1,2$.
We show the following:
\begin{itemize}
\item[(i)] For the case when $(k_1-1)(k_2-1) >0$, we take $k_1=k_2=\infty$ and show that the solution $u$ to \eqnref{NCE} with $f$ defined above satisfies
\begin{equation}\label{reverse1}
|\nabla u(z)| \gtrsim \epsilon^{-1/2}
\end{equation}
for some $z \in \mathbb{R}^2 \setminus \overline{D}$.
\item[(ii)] For the case when $(k_1-1)(k_2-1) < 0$, we take either $k_1=0, \ k_2=\infty$ or $k_1=\infty, \ k_2=0$ and show that the solution $u$ to \eqnref{NCE} with $f$ defined above satisfies
\begin{equation}\label{reverse2}
|\nabla^2 u(z)| \gtrsim \epsilon^{-1/2}
\end{equation}
for some $z \in \mathbb{R}^2 \setminus \overline{D}$. ($\nabla u$ is bounded.)
\end{itemize}
Similar estimates hold for the solution to the homogeneous problem with $H(x)=x_1$.

According to the definition \eqnref{hjdef} of $h_j$, we have
\begin{equation}\label{hjexample}
h_j(\zeta)= \frac{\beta}{\zeta-1}, \quad \zeta \in D_j^*.
\end{equation}

Suppose that $k_1=k_2=\infty$. Then, $\lambda_1+\lambda_2=1$ and $\lambda_1-\lambda_2=0$.
Since $4 \lambda_1\lambda_2=1$ in this case, it follows from the formulas \eqnref{A1def} and \eqnref{A2def} of $A_j$ that
$$
V(\zeta) = A_1(\zeta) + A_2(\zeta) = w_1(R_1^2 \bar{\Gz}^{-1}) - w_1(\rho^2 \zeta) + w_2(R_2^2 \zeta^{-1}) - w_2(\rho^{-2} \bar{\Gz})
$$
if $R_1 < |\zeta| \leq R_2$. Then the formulas \eqnref{w1def2} and \eqnref{w2def2} for $w_j$ yield
\begin{align*}
\partial_\zeta V(\zeta) &= - \rho^2 w_1'(\rho^2 \zeta) - R_2^2 \zeta^{-2} w_2'(R_2^2 \zeta^{-1}) \\
& = -\sum_{l=0}^\infty \rho^{2l+2} h_1' (\rho^{2l+2}\zeta) - \sum_{l=0}^\infty \rho^{-2l}R_2^2 \zeta^{-2} h_2' (\rho^{-2l}R_2^2 \zeta^{-1}) .
\end{align*}
It thus follows from \eqnref{hjexample} that
\begin{equation}\label{pzV}
\partial_\zeta V(\zeta) = \beta \sum_{l=0}^\infty \left[ \frac{\rho^{2l+2}}{(1-\rho^{2l+2}\zeta)^2} + \frac{\rho^{2l} R_2^{-2}}{(1-\rho^{2l} R_2^{-2}\zeta)^2} \right] .
\end{equation}
In particular, we have
$$
\partial_\zeta V(-1) = \beta \sum_{l=0}^\infty \left[ \frac{\rho^{2l+2}}{(1+\rho^{2l+2})^2} + \frac{\rho^{2l} R_2^{-2}}{(1+\rho^{2l} R_2^{-2})^2} \right] \gtrsim \frac{1}{1-\rho} \gtrsim \frac{1}{\beta}.
$$
Note that
\begin{equation}\label{dGzdz}
\frac{\partial^k \zeta}{\partial z^k}= (-1)^k \frac{k! \beta}{z^{k+1}} = (-1)^k k! \frac{(\zeta-1)^{k+1}}{\beta^{k}}, \quad k=1,2,\ldots.
\end{equation}
Thus we have
\begin{equation}
\big| \partial_z (V\circ T)(T^{-1}(-1)) \big| = \Big| \partial_\zeta V(-1) \frac{\partial \zeta}{\partial z} \Big|_{\zeta=-1} \Big| \gtrsim \frac{1}{\beta},
\end{equation}
and we arrive at \eqnref{reverse1}.

\medskip

We now consider the case when $(k_1-1)(k_2-1) < 0$. If $k_1=0$ and $k_2=\infty$, then $\lambda_1=-1/2$ and $\lambda_2=1/2$, and hence $\lambda_1+\lambda_2=0$, $\lambda_1-\lambda_2=-1$, and $4 \lambda_1\lambda_2=-1$. We have from \eqnref{A1def} and \eqnref{A2def} that
$$
V(\zeta) = -w_1(\rho\zeta) - w_1(\rho^2 \zeta) + w_2(\rho^{-1} \bar{\Gz}) - w_2(\rho^{-2} \bar{\Gz})
$$
if $R_1 < |\zeta| \leq R_2$. Thus we have
\begin{equation}
\partial_\zeta V(\zeta) = - \rho w_1'(\rho\zeta) - \rho^2 w_1'(\rho^2 \zeta).
\end{equation}

We see from \eqnref{w1def2} and \eqnref{hjexample} that
$$
w_1(\zeta) = \beta \sum_{l=0}^\infty (-1)^{l+1} \frac{1}{1-\rho^{2l}\zeta} = -\beta \sum_{k=0}^\infty \frac{\zeta^k}{1+\rho^{2k}},
$$
and hence
$$
w_1'(\zeta) = -\beta \sum_{k=1}^\infty \frac{k\zeta^{k-1}}{1+\rho^{2k}}.
$$
Thus, we have
$$
\partial_\zeta V(\zeta) = \beta \sum_{k=1}^\infty \frac{k (\rho^{k} + \rho^{2k}) \zeta^{k-1}}{1+\rho^{2k}} .
$$
We also have
$$
\partial_\zeta^2 V(\zeta) = \beta \sum_{k=2}^\infty \frac{k(k-1) (\rho^{k} + \rho^{2k}) \zeta^{k-2}}{1+\rho^{2k}} .
$$
Note that
$$
\partial_z^2 (V\circ T)(z) = \partial_\zeta^2 V(\zeta) \left(\frac{\partial \zeta}{\partial z} \right)^2 + \partial_\zeta V(\zeta) \frac{\partial^2 \zeta}{\partial z^2}.
$$
Thus \eqnref{dGzdz} yields
$$
\partial_z^2 (V\circ T)(z) = \frac{(\zeta-1)^3}{\beta} \Big[ \sum_{k=1}^\infty \frac{k(k+1) (\rho^{k} + \rho^{2k}) \zeta^{k-1}}{1+\rho^{2k}} +  \sum_{k=2}^\infty \frac{k(k-1) (\rho^{k} + \rho^{2k}) \zeta^{k-2}}{1+\rho^{2k}} \Big].
$$
Note that
$$
\int_{-\pi}^{\pi} \beta \partial_z^2 (V\circ T)(T^{-1}(e^{i\theta})) d\theta= - 2\pi C,
$$
where
$$
C= \frac{2(\rho + \rho^{2})}{1+\rho^{2}} +  \frac{2(\rho^{2} + \rho^{4})}{1+\rho^{4}}.
$$
It means that there is a point, say $\zeta_*$, such that $|\zeta_*|=1$ and
$$
\beta \partial_z^2 (V\circ T)(T^{-1}(\zeta_*)) \le - C.
$$
Note that $C \ge 1$ if $\epsilon$ is sufficiently small. Thus we have
\begin{equation}\label{7100}
\big| \partial_z^2 (V\circ T)(T^{-1}(\zeta_*)) \big| \gtrsim \frac{1}{\beta} ,
\end{equation}
which yields \eqnref{reverse2}.

If $k_1=\infty$ and $k_2=0$, then $\lambda_1+\lambda_2=0$, $\lambda_1-\lambda_2=1$, and $4 \lambda_1\lambda_2=-1$. We have from \eqnref{A1def} and \eqnref{A2def} that
$$
V(\zeta) = w_1(\rho\zeta) - w_1(\rho^2 \zeta) - w_2(\rho^{-1} \bar{\Gz}) - w_2(\rho^{-2} \bar{\Gz})
$$
if $R_1 < |\zeta| \leq R_2$. Thus we have
\begin{equation}
\partial_{\bar{\Gz}} V(\zeta) = - \rho^{-1} w_2'(\rho^{-1} \bar{\Gz}) - \rho^{-2} w_2'(\rho^{-2} \bar{\Gz}).
\end{equation}
Now in the same way as above, one can show that there is a $\zeta_{**}$ such that $|\zeta_{**}|=1$ and
\begin{equation}\label{7200}
\big| \partial_{\bar{z}}^2 (V\circ T)(T^{-1}(\zeta_{**})) \big| \gtrsim \frac{1}{\beta} ,
\end{equation}
which again yields \eqnref{reverse2}.