# Interaction energies in paranematic colloids

Dmitry Golovaty Department of Mathematics, University of Akron.
<EMAIL_ADDRESS>Jamie Taylor Departamento de Métodos Cuantitativos, CUNEF
Universidad<EMAIL_ADDRESS>Raghavendra Venkatraman Courant Institute
of Mathematical Sciences, New York University<EMAIL_ADDRESS>Arghir
Zarnescu Basque Center for Applied Mathematics, Ikerbasque Foundation and
”Simion Stoilow” Institute of the Romanian Academy<EMAIL_ADDRESS>
###### Abstract

We consider a system of colloidal particles embedded in a paranematic—an
isotropic phase of a nematogenic medium above the temperature of the nematic-
to-isotropic transition. In this state, the nematic order is induced by the
boundary conditions in a narrow band around each particle and it decays
exponentially in the bulk.

We develop rigorous asymptotics of the linearization of the appropriate
variational model that allow us to describe weak far-field interactions
between the colloidal particles in two dimensional paranematic suspensions. We
demonstrate analytically that decay rates of solutions to the full nonlinear
and linear problems are similar and verify numerically that the interactions
between the particles in these problems have similar dependence on the
distance between the particles. Finally, we perform Monte-Carlo simulations
for a system of colloidal particles in a paranematic and describe the
statistical properties of this system.

###### Contents

  1. 1 Introduction and main results
  2. 2 Two particles and constant boundary conditions
    1. 2.1 Interaction energies for large vanishing separation between particles
    2. 2.2 Interaction energies for O(1) separation between particles
  3. 3 Two particles and nonconstant boundary conditions
  4. 4 Interaction energies of multiple particles
  5. 5 Numerics and comparison to nonlinear models
  6. 6 Monte Carlo simulation of the pairwise energy
  7. 7 Acknowledgements
  8. A Sobolev spaces and trace theory
  9. B The case of a single particle
  10. C Estimates of modified Bessel functions of the second kind
  11. D Some facts about Chebyshev polynomials of the first and second kinds
  12. E  A brief introduction to the Landau-de Gennes model

## 1\. Introduction and main results

We aim to initiate the study of interaction energies between colloidal
particles in a nematic liquid crystal environment. There exists a significant
body of physics literature on this topic (see for instance [14, 13, 15]),
mostly based on simulations of certain variational models. The analytical
intuition behind these interactions follows ideas developed in the seminal
paper [12]. In this work, interactions between colloidal particles are
established based on a suitable linearisation at infinity and on formal
analogies with a classical theme, namely interactions of electrostatic
multipoles. A rigorous understanding of these interactions is still missing in
the case of several particles, while for the case of a single particle was
considered in the recent work [2].

The main goal in this paper is to provide a rigorous underpinning to the
intuition developed in the physical literature, aiming to obtain explicit
estimates quantifying the interaction in the case of several particles
expressed in terms of the geometric and material parameters of the problem.
Models typically used describe this physical setting are nonlinear but
following formal ideas in [4, 16], we reduce the problem to the linearization
around the isotropic state and discuss the precise analytical meaning in which
the solution to the resulting linear problem approximate the minimizers of the
nonlinear problem. Further, we conduct the detailed analytical study of the
linear problem, and then show by numerical experiments that the nonlinear
version of the problem shares a number of qualitative features of our linear
analysis.

In the long-term we will be interested in understanding a Landau-de Gennes
model of nematic liquid crystals. The main modelling features of this model
based on a tensor-valued order parameter are presented in Appendix E. In the
current paper we focus on a simpler, vector-valued model that retains the
relevant features of the Landau-de Gennes approach. To this end, suppose that
the liquid crystal is described by $u:\Omega\to\mathbb{R}^{2}$, where
$\Omega\subset\mathbb{R}^{2}$ is an open, smooth and not necessarily bounded
domain that models the container containing the liquid crystal material. Let

$W(u)=k(T){|u|}^{2}-2{|u|}^{4}+{|u|}^{6}$ (1.1)

be a bulk potential, the minima of which describe some physical system that
may undergo a phase transition at some critical temperature $T=T^{*}$, i.e. in
mathematical terms, the type and number of minima change at this temperature.

An examination of $W$ reveals that it has exactly one minimum at the isotropic
state $u=0$ when $k(T)>\frac{4}{3}.$ When $\frac{4}{3}\geq k(T)>1$ there is a
global minimum at $u=0$ and a local minimum at $|u|=\alpha(T)$ that represents
a metastable ordered state. When $k(T)=1$ both minima have equal depth, while
$|u|=\alpha(T)$ and $u=0$ become the global minimum and a local minimum,
respectively, when $1>k(T)>0.$ When $k(T)<0,$ the circle $|u|=\alpha(T)$ is
the global minimal set, while $u=0$ is the local maximum of $W$. We call the
temperature $T^{*}$ satisfying $k(T^{*})=1$ the temperature of the phase
transition between the isotropic and the ordered states.

We study the minimizers of the functional

$E_{\varepsilon}(u)=\int_{\Omega}{\frac{1}{2}{|\nabla
u|}^{2}+\frac{1}{4\varepsilon^{4}}W(u)},$ (1.2)

where we assume that $u\in H^{1}(\Omega)$ satisfies Dirichlet boundary data on
$\partial\Omega$. We will be focused on domains $\Omega$ that are exterior to
a collection of colloidal particles, and seek to understand their inter-
particle interactions as mediated by the background ordered state. As a first
step in this program, we are interested in a so-called _paranematic_ regime
when $k(T)>\frac{4}{3}.$ In this regime, the potential $W$ is convex and has a
single minimum at the isotropic state $u=0.$ Nonetheless, whilst the ground
state is isotropic in the bulk, some residual nematic ordering may be induced
by the boundary conditions.

Note that interactions between spherical particles immersed in an isotropic
phase of a nematogenic fluid were investigated in the physical literature in
the past by considering formal asymptotics [4]-[16] in three dimensions. Here
we will focus instead on rigorous understanding of the regime of two spherical
colloids in $\mathbb{R}^{2}$ when the domain is the whole space. To fix ideas,
we suppose thus that there are two identical, spherical colloidal particles
$\widehat{B}_{1}^{\varepsilon}$ and $\widehat{B}^{\varepsilon}_{2},$ each of
radius $1$, that are separated by distance $2b\varepsilon^{2}>0.$ We also take
$\Omega=\mathbb{R}^{2}\setminus\overline{\widehat{B}_{1}^{\varepsilon}\cup\widehat{B}_{2}^{\varepsilon}}$.

We show in Proposition 2.6 that in the paranematic regime, the unique solution
to the nonlinear Euler-Lagrange problem

$\displaystyle\Delta u$
$\displaystyle=\frac{1}{4\varepsilon^{4}}\nabla_{u}W(u),\quad u\in
H^{1}\mathopen{}\mathclose{{}\left(\mathbb{R}^{2}\setminus\overline{\widehat{B}_{1}^{\varepsilon}\cup\widehat{B}_{2}^{\varepsilon}}}\right)$
(1.3) $\displaystyle u$ $\displaystyle=g_{1}\qquad\mbox{ on
}\partial\widehat{B}_{1}^{\varepsilon}$ $\displaystyle u$
$\displaystyle=g_{2}\qquad\mbox{ on }\partial\widehat{B}_{2}^{\varepsilon}\,.$

for (1.2) has the same rate of decay as that of a solution of the
corresponding linearization of (1.3) around the state $u\equiv 0,$ namely

$\displaystyle\Delta u$ $\displaystyle=\frac{k(T)}{2\varepsilon^{4}}u,\quad
u\in
H^{1}\mathopen{}\mathclose{{}\left(\mathbb{R}^{2}\setminus\overline{\widehat{B}_{1}^{\varepsilon}\cup\widehat{B}_{2}^{\varepsilon}}}\right)$
(1.4) $\displaystyle u$ $\displaystyle=g_{1}\qquad\mbox{ on
}\partial\widehat{B}_{1}^{\varepsilon}$ $\displaystyle u$
$\displaystyle=g_{2}\qquad\mbox{ on }\partial\widehat{B}_{2}^{\varepsilon}\,.$

This observation allows us to conjecture that far-field paranematic-mediated
interactions between two particles in a nematogenic medium should depend on
the distance between the particles in a way similar to that for the particles
in the corresponding linear problem.

By rescaling $\varepsilon$, we can assume that $k(T)=2$ in (1.1), hence the
linear PDE we will consider in the sequel is

$\displaystyle\Delta U$ $\displaystyle=\frac{1}{\varepsilon^{4}}U,\quad U\in
H^{1}\mathopen{}\mathclose{{}\left(\mathbb{R}^{2}\setminus\overline{\widehat{B}_{1}^{\varepsilon}\cup\widehat{B}_{2}^{\varepsilon}}}\right)$
(1.5) $\displaystyle U$ $\displaystyle=G_{1}\qquad\mbox{ on
}\partial\widehat{B}_{1}^{\varepsilon}$ $\displaystyle U$
$\displaystyle=G_{2}\qquad\mbox{ on }\partial\widehat{B}_{2}^{\varepsilon}\,.$

It is clear that for the linear PDE (1.5), the vectorial nature of $U$ is
unimportant and, therefore, whenever possible we will assume $U$ to be a
scalar.

For concreteness, we set $\widehat{B}_{1}^{\varepsilon}$ to be an open unit
ball centered at $(0,1+b\varepsilon^{2}),$ and $\widehat{B}_{2}^{\varepsilon}$
an open unit ball centered at $(0,-1-b\varepsilon^{2}),$ so that in this
scaling the distance of separation between the balls is $2b\varepsilon^{2}>0.$
This choice of geometry is motivated by the desire to compute the energy of
interaction between the balls mediated by their paranematic surroundings.

It will often be convenient to rescale the problem via the change of variables
$x^{\prime}=\frac{x}{\varepsilon^{2}},y^{\prime}=\frac{y}{\varepsilon^{2}}$
and subsequently dropping the primes. Then, setting $B_{1}^{\varepsilon}$ to
be the ball of radius $\frac{1}{\varepsilon^{2}}$ centered at
$(0,\frac{1}{\varepsilon^{2}}+b)$ and $B_{2}^{\varepsilon}$ to be the ball of
radius $\frac{1}{\varepsilon^{2}}$ centered at
$(0,-\frac{1}{\varepsilon^{2}}-b),$ we have

$\displaystyle\Delta u$ $\displaystyle=u,\quad\quad u\in
H^{1}\mathopen{}\mathclose{{}\left(\mathbb{R}^{2}\setminus\overline{B_{1}^{\varepsilon}\cup
B_{2}^{\varepsilon}}}\right)$ (1.6) $\displaystyle u$
$\displaystyle=g_{1}\qquad\mbox{ on }\partial B_{1}^{\varepsilon}$
$\displaystyle u$ $\displaystyle=g_{2}\qquad\mbox{ on }\partial
B_{2}^{\varepsilon}\,.$

In the above, functions with lower case letters $(u,g_{1},g_{2})$ represent
scaled versions of their upper-case counterparts $(U,G_{1},G_{2})$. We point
out that in the blown up variables, the separation between the balls is
$2b>0.$

For the rescaling as above, we observe that the natural quadratic energies
associated to the two settings are equal:

$\int_{\mathbb{R}^{2}\setminus\overline{\widehat{B}_{1}^{\varepsilon}\cup\widehat{B}_{2}^{\varepsilon}}}\biggl{(}|\nabla
U|^{2}+\frac{1}{\varepsilon^{4}}U^{2}\biggr{)}\,dx=\int_{\mathbb{R}^{2}\setminus\overline{B_{1}^{\varepsilon}\cup
B_{2}^{\varepsilon}}}\biggl{(}|\nabla u|^{2}+u^{2}\biggr{)}\,dx\,.$

Note that, although formally these energies can be thought of as a leading
order approximation of (1.2) when the supremum norm of $u$ is small, this is
not true in the current case as the boundary data $g_{i},$ $i=1,2$ is of order
$1.$

It will be shown by direct energy comparison with a competitor, that the
energy of the unique solution to (1.6) satisfies

$\int_{\mathbb{R}^{2}\setminus\overline{B_{1}^{\varepsilon}\cup
B_{2}^{\varepsilon}}}\biggl{(}|\nabla
u|^{2}+u^{2}\biggr{)}\,dx\leqslant\frac{C}{\varepsilon^{2}}\|g\|_{L^{2}}^{2}\,.$
(1.7)

For the linear problem, the goal in this paper is to give a precise energy
expansion of the first two terms of the minimum energy $E(u),$ in terms of the
parameter $\varepsilon,$ and quantify the expansion. To be precise, let us
note that the problems (1.5) and (1.6) are associated with variational
principles, and the solutions to these PDEs arise as unique minimizers of
strictly convex energies. Focusing on (1.6), the associated energy is

$F_{\varepsilon}(u):=\frac{1}{2}\int_{\Omega_{\varepsilon}}\Bigl{(}|\nabla
u|^{2}+u^{2}\Bigr{)}\,dx\,,\quad\Omega_{\varepsilon}:=\mathbb{R}^{2}\setminus(B_{1}^{\varepsilon}\cup
B_{2}^{\varepsilon})\,.$ (1.8)

It is easy to see that since we are working in two dimensions, the energy
associated with (1.5) is equal to that with (1.6). We set

$\kappa_{\varepsilon}(g_{1},g_{2}):=\min\Bigl{\\{}F_{\varepsilon}(w):w\in
H^{1}(\Omega_{\varepsilon}),w=g_{i}\,\,\mbox{ on }\partial
B_{i}^{\varepsilon}\Bigr{\\}}\,.$ (1.9)

Our first main result, to be provided in Section 2 concerns the first two
terms in an asymptotic expansion of $\kappa_{\varepsilon}(g_{1},g_{2})$ in
powers of $\varepsilon.$ and for constant boundary conditions
$g_{1},g_{2}\in\mathbb{R}$. The leading terms are of order
$O(\varepsilon^{-2})$, and correspond to the energy of each individual
particle. However, due to the presence of two particles (rather than one), the
particles interact, and there is a correction to the leading order energy
which occurs at order $O(\varepsilon^{-1}).$ Exactly computing this
interaction energy is the main contribution of our work. More precisely, in
Subsection 2 we will provide:

######  Theorem 1.1.

Let
$B_{1}^{\varepsilon}=\\{x\in\mathbb{R}^{2},d(x,(0,\varepsilon^{-2}+b))<\varepsilon^{-2}\\}$
and
$B_{2}^{\varepsilon}=\\{x\in\mathbb{R}^{2},d(x,(0,-\varepsilon^{-2}-b))<\varepsilon^{-2}\\}$
for $\varepsilon\in(0,\frac{1}{2}]$ and $b\gg 1.$

For constant boundary conditions $g_{1},g_{2}\in\mathbb{R}$ we consider the
energy $F_{\varepsilon}$ defined as in (1.8) and its minimum $k_{\varepsilon}$
as defined in (1.9). Then we have the estimate

$\mathopen{}\mathclose{{}\left|\kappa_{\varepsilon}-\frac{2\pi}{\varepsilon^{2}}(g_{1}^{2}+g_{2}^{2})+4g_{1}g_{2}\frac{e^{-2b}\sqrt{\pi}}{\varepsilon}}\right|\leqslant\frac{Ce^{-4b}}{\varepsilon}.$
(1.10)

######  Remark 1.2.

In the other regime we will consider, in the limit $\varepsilon\to 0^{+}$ for
fixed separation $b=O(1)$ as $\varepsilon\to 0^{+}$, we will only provide the
formal asymptotics, in Section 2.2 , in order to limit the size of the
presentation, although rigorous statements should be obtainable in this case
as well. More specifically in this case we have

$\kappa_{\varepsilon}=\frac{2\pi}{\varepsilon^{2}}\mathopen{}\mathclose{{}\left(g_{1}^{2}+g_{2}^{2}}\right)+\frac{\sqrt{\pi/2}}{\varepsilon}\mathopen{}\mathclose{{}\left[2e^{-4b}+\frac{1}{\sqrt{2}}e^{-8b}}\right](g_{1}^{2}+g_{2}^{2})\\\
-\frac{4\sqrt{\pi}}{\varepsilon}\mathopen{}\mathclose{{}\left[e^{-2b}+\frac{1}{\sqrt{3}}e^{-6b}}\right]g_{1}g_{2}+O(1),$

Moving on to the more general case of non-constant boundary conditions, we
will show in Section 3 the following:

######  Theorem 1.3.

Let
$B_{1}^{\varepsilon}=\\{x\in\mathbb{R}^{2},d(x,(0,\varepsilon^{-2}+b))<\varepsilon^{-2}\\}$
and
$B_{2}^{\varepsilon}=\\{x\in\mathbb{R}^{2},d(x,(0,-\varepsilon^{-2}-b))<\varepsilon^{-2}\\}$
for $\varepsilon\in(0,\frac{1}{2}]$ and $b\gg 1.$

For $g_{1},g_{2}$ smooth functions on $\partial B_{1}^{\varepsilon}$
respectively $\partial B_{2}^{\varepsilon}$, we consider the energy
$F_{\varepsilon}$ defined as in (1.8) and its minimum $k_{\varepsilon}$ as
defined in (1.9). Then we have the estimate

$\displaystyle\mathopen{}\mathclose{{}\left|\kappa_{\varepsilon}-\frac{2\pi}{\varepsilon^{2}}\Biggl{(}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{\partial
B_{1}}|g_{1}|^{2}\,d\mathcal{H}^{1}+\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{\partial
B_{2}}|g_{2}|^{2}\,d\mathcal{H}^{1}\Biggr{)}+2\frac{e^{-2b}\sqrt{\pi}}{\varepsilon}\Bigl{(}g_{1}(p)\overline{g_{2}(q)}+\overline{g_{1}(p)}g_{2}(q)\Bigr{)}}\right|$
$\displaystyle\quad\quad\leqslant\frac{Ce^{-4b}(\|g_{1}\|_{H^{1}(\partial
B_{1})}^{2}+\|g_{2}\|_{H^{1}(\partial B_{2})}^{2})}{\varepsilon}\,,$

where $p\in\partial B_{1}$ and $q\in\partial B_{2}$ denote the points on the
two respective circles that are closest to each other.

Furthermore, we will briefly explain in Section 4 how the results for two-
particle interactions can be extended to the case of several particles, under
suitable assumptions.

In Section 5 through numerics we will explore the similarities (and
differences) at the level of minimizers and energy scaling, when the simple
quadratic potential $\frac{1}{2}u^{2}$ in the energy $F_{\varepsilon}(u)$ is
replaced by the nonlinear (and the more physical) potential $W(u)$
corresponding to the paranematic setting.

The results presented here are providing analytical first steps towards a
rigorous understanding of multiple particles interactions and open the doors
of a fascinating world. In order to offer a glimpse of the future explorations
we provide in Figure 1 some Monte Carlo simulations based on the ideas
developed here. These show configurations of several particles with boundary
conditions having different topological degrees and the details are provided
in Section 6.

(a) Degree $2$

(b) Degree $3$

(c) Degree $5$

Figure 1: Some Monte-Carlo simulations for multiple particles

## 2\. Two particles and constant boundary conditions

The focus of this section is on the case when the boundary conditions in (1.6)
are constant, so that $g_{1},g_{2}\in\mathbb{R}.$ Recall that
$\Omega_{\varepsilon}=\mathbb{R}^{2}\setminus\overline{B_{1}^{\varepsilon}\cup
B_{2}^{\varepsilon}}$ is the domain exterior to two large balls of the radius
$1/\varepsilon^{2}$ each and situated at distance $2b$ away from each other.
We will provide an energy expansion in two cases:

  1. 1.

that holds for a _fixed_ $\varepsilon\in(0,\frac{1}{2}]$ and $b,$ in the limit
of $b\gg 1,$ and

  2. 2.

that holds in the limit $\varepsilon\to 0^{+}$ for fixed separation $b=O(1)$
as $\varepsilon\to 0^{+}.$

The first case, to be treated in the next subsection will be studied
rigorously, providing all the details, while for the other case, to be treated
in Subsection 2.2 we will only provide the formal asymptotics, in order to
limit to size of the presentation.

###  Interaction energies for large vanishing separation between particles

The main result of this section is an expansion of the energy of $U$ in terms
of $\varepsilon$, namely the one stated in Theorem 1.1.

The proof of Theorem 1.1 is contained in a sequence of Lemmas. For $i=1,2,$
let

$\Psi_{i}(x):=\frac{K_{0}\big{(}|x-a_{i}^{\varepsilon}|\big{)}}{K_{0}\big{(}\frac{1}{\varepsilon^{2}}\big{)}}\,,\quad
x\in\Omega_{\varepsilon}$ (2.1)

where $K_{0}$ is a modified Bessel function of the second kind (see Appendix C
for details). One can check that this is the solution of the exterior problem

$\displaystyle\Delta\Psi_{i}=\Psi_{i}\quad\quad\mbox{ in
}\mathbb{R}^{2}\backslash B_{i}^{\varepsilon},$ (2.2)
$\displaystyle\Psi_{i}=1\quad\quad\quad\ \mbox{ on }\partial
B_{i}^{\varepsilon}$

Let us define $\sigma:\\{1,2\\}\to\\{1,2\\}$ via $\sigma(1)=2,\sigma(2)=1.$
Suppose that $Z_{i}\in
H^{1}\mathopen{}\mathclose{{}\left(\Omega_{\varepsilon}}\right)$ solve

$\displaystyle\Delta Z_{i}=Z_{i}\quad\quad\mbox{ in }\Omega_{\varepsilon},$
(2.3) $\displaystyle Z_{i}=0\quad\ \quad\quad\mbox{ on }\partial
B_{i}^{\varepsilon},$ $\displaystyle Z_{i}=-\Psi_{i}\quad\quad\mbox{ on
}\partial B_{\sigma(i)}^{\varepsilon}$

for $i=1,2.$ Then it is easily seen that the unique solution to (1.6) is given
by

$\displaystyle U=g_{1}(\Psi_{1}+Z_{1})+g_{2}(\Psi_{2}+Z_{2}).$ (2.4)

If for any $\xi_{1},\xi_{2}\in H^{1}(\Omega_{\varepsilon})$ we define

$\displaystyle\langle\xi_{1},\xi_{2}\rangle_{\varepsilon}:=\langle\xi_{1},\xi_{2}\rangle_{H^{1}(\Omega_{\varepsilon})}=\int_{\Omega_{\varepsilon}}\xi_{1}\xi_{2}+\nabla\xi_{1}\cdot\nabla\xi_{2}$

and denote $\|\xi\|_{\varepsilon}^{2}:=\langle\xi,\xi\rangle_{\varepsilon},$
we observe that

$\displaystyle F_{\varepsilon}(U)=\frac{1}{2}\|U\|_{\varepsilon}^{2}.$

Then,

$\displaystyle
2F_{\varepsilon}(U)=\|U\|_{\varepsilon}^{2}=g_{1}^{2}\|\Psi_{1}+Z_{1}\|_{\varepsilon}^{2}+2g_{1}g_{2}\langle\Psi_{1}+Z_{1},\Psi_{2}+Z_{2}\rangle_{\varepsilon}+g_{2}^{2}\|\Psi_{2}+Z_{2}\|_{\varepsilon}^{2}$
(2.5)
$\displaystyle\qquad=\begin{pmatrix}g_{1}&g_{2}\end{pmatrix}\begin{pmatrix}\|\Psi_{1}+Z_{1}\|_{\varepsilon}^{2}&\langle\Psi_{1}+Z_{1},\Psi_{2}+Z_{2}\rangle_{\varepsilon}\\\
\langle\Psi_{1}+Z_{1},\Psi_{2}+Z_{2}\rangle_{\varepsilon}&\|\Psi_{2}+Z_{2}\|_{\varepsilon}^{2}\end{pmatrix}\begin{pmatrix}g_{1}\\\
g_{2}\end{pmatrix}.$

Our first lemma expresses each of the terms in the above matrix in terms of
certain boundary integrals. Naturally, this is done using integration by
parts– for this purpose we let $\nu_{i},i=1,2$ denote unit normals that point
_towards the centers of the discs_ $B_{i}^{\varepsilon}.$ In particular, for
the exterior domain $\Omega_{\varepsilon},$ these represent outward unit
normals. We will collectively refer to both these normals (i.e., as outward
unit normal to $\partial\Omega_{\varepsilon}$) via $\nu$.

######  Lemma 2.1.

We have the following identities

$\displaystyle\|\Psi_{2}+Z_{2}\|_{\varepsilon}^{2}=\|\Psi_{1}+Z_{1}\|_{\varepsilon}^{2}=-\frac{2\pi}{\varepsilon^{2}}\frac{K_{0}^{\prime}\bigl{(}\frac{1}{\varepsilon^{2}}\bigr{)}}{K_{0}\bigl{(}\frac{1}{\varepsilon^{2}}\bigr{)}}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\biggl{(}\frac{\partial\Psi_{1}}{\partial\nu_{2}}+\frac{\partial
Z_{1}}{\partial\nu_{2}}\biggr{)}\,d\mathcal{H}^{1},$ (2.6)

and

$\displaystyle\langle\Psi_{1}+Z_{1},\Psi_{2}+Z_{2}\rangle_{\varepsilon}=\int_{\partial
B_{2}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{2}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}.$ (2.7)

###### Proof.

To prove (2.6), we compute $\|\Psi_{1}+Z_{1}\|_{\varepsilon}^{2},$ with the
other term being symmetrical. Computing, and using the PDE and boundary
conditions satisfied by $\Psi_{1}$ and $Z_{1},$ we observe

$\displaystyle\langle\Psi_{1}+Z_{1}\rangle_{\varepsilon}^{2}=\|\Psi_{1}\|_{\varepsilon}^{2}+2\langle\Psi_{1},Z_{1}\rangle_{\varepsilon}+\|Z_{1}\|_{\varepsilon}^{2}$
(2.8)
$\displaystyle\quad\quad=\int_{\Omega_{\varepsilon}}\mathopen{}\mathclose{{}\left(\Psi_{1}^{2}+|\nabla\Psi_{1}|^{2}+2\Psi_{1}Z_{1}+2\nabla\Psi_{1}\cdot\nabla
Z_{1}+|Z_{1}|^{2}+|\nabla Z_{1}|^{2}}\right)\,dx$
$\displaystyle\quad\quad=\int_{\Omega_{\varepsilon}}\Psi_{1}\mathopen{}\mathclose{{}\left(\Psi_{1}-\Delta\Psi_{1}}\right)\,dx+\int_{\partial\Omega_{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{1}}{\partial\nu}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad\
+2\int_{\Omega_{\varepsilon}}Z_{1}\mathopen{}\mathclose{{}\left(\Psi_{1}-\Delta\Psi_{1}}\right)\,dx+2\int_{\partial\Omega_{\varepsilon}}Z_{1}\frac{\partial\Psi_{1}}{\partial\nu}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad\ +\int_{\Omega_{\varepsilon}}Z_{1}(Z_{1}-\Delta
Z_{1})\,dx+\int_{\partial\Omega_{\varepsilon}}Z_{1}\frac{\partial
Z_{1}}{\partial\nu}\,d\mathcal{H}^{1}$ $\displaystyle\quad\quad=\int_{\partial
B_{1}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}+\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-2\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial
Z_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad=\int_{\partial
B_{1}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\mathopen{}\mathclose{{}\left(\frac{\partial\Psi_{1}}{\partial\nu_{2}}+\frac{\partial
Z_{1}}{\partial\nu_{2}}}\right)\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad=-\frac{2\pi}{\varepsilon^{2}}\frac{K_{0}^{\prime}\mathopen{}\mathclose{{}\left(\frac{1}{\varepsilon^{2}}}\right)}{K_{0}\mathopen{}\mathclose{{}\left(\frac{1}{\varepsilon^{2}}}\right)}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\mathopen{}\mathclose{{}\left(\frac{\partial\Psi_{1}}{\partial\nu_{2}}+\frac{\partial
Z_{1}}{\partial\nu_{2}}}\right)\,d\mathcal{H}^{1}.$

Similarly, to prove (2.7), we notice

$\displaystyle\langle\Psi_{1}+Z_{1},\Psi_{2}+Z_{2}\rangle_{\varepsilon}=\int_{\Omega_{\varepsilon}}(\Psi_{1}+Z_{1})(\Psi_{2}+Z_{2})+(\nabla\Psi_{1}+\nabla
Z_{1})\cdot(\nabla\Psi_{2}+\nabla Z_{2})\,dx$ (2.9)
$\displaystyle\quad=\int_{\Omega_{\varepsilon}}(\Psi_{1}-\Delta\Psi_{1})\Psi_{2}\,dx+\int_{\partial\Omega_{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu}\Psi_{2}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad+\int_{\Omega_{\varepsilon}}(\Psi_{1}-\Delta\Psi_{1})Z_{2}\,dx+\int_{\partial\Omega_{\varepsilon}}Z_{2}\frac{\partial\Psi_{1}}{\partial\nu}\,\mathcal{H}^{1}$
$\displaystyle\quad\quad+\int_{\Omega_{\varepsilon}}(\Psi_{2}-\Delta\Psi_{2})Z_{1}\,dx+\int_{\partial\Omega_{\varepsilon}}Z_{1}\frac{\partial\Psi_{2}}{\partial\nu}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad+\int_{\Omega_{\varepsilon}}(Z_{1}-\Delta
Z_{1})Z_{2}\,dx+\int_{\partial\Omega_{\varepsilon}}Z_{2}\frac{\partial
Z_{1}}{\partial\nu}\,d\mathcal{H}^{1}$ $\displaystyle\quad=\int_{\partial
B_{1}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{1}}\Psi_{2}\,d\mathcal{H}^{1}+\int_{\partial
B_{2}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad-\int_{\partial
B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial\Psi_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{2}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}$ $\displaystyle\quad=\int_{\partial
B_{2}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{2}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}.$

∎

Next, we have a lemma that controls the normal derivative of the function $Z$
in $H^{-1/2}(\partial\Omega_{\varepsilon})$ by the energy. The underlying
subtlety, is of course, that the domain $\Omega_{\varepsilon}$ varies in
$\varepsilon,$ and we must obtain estimates that are uniform in $\varepsilon.$

######  Lemma 2.2.

The functions $Z_{j}$ have $H^{1}$ norms bounded by

$\displaystyle\biggl{\\{}\int_{\Omega_{\varepsilon}}|Z_{j}|^{2}+|\nabla
Z_{j}|^{2}\,dx\biggr{\\}}^{\nicefrac{{1}}{{2}}}\leqslant\frac{C}{\sqrt{\varepsilon}}K_{0}(2b),$
(2.10)

where $C$ is independent of $\varepsilon.$

###### Proof.

We first make the observation that the functions $Z_{j}$ satisfying the PDE
(2.3) are the unique minimizers of the $H^{1}$ norm, subject to their own
boundary conditions. Therefore, the desired estimate follows by the
construction of a competitor and comparing energies. Without loss of
generality, we fix $j=1$. Our competitor $\zeta\in
H^{1}(\Omega_{\varepsilon})$ must be constructed satisfying the boundary
conditions for $\zeta=Z_{1}$ on $\partial\Omega_{\varepsilon},$ so that
$\zeta=0$ on $\partial B_{1}$ and $\zeta=-\Psi_{1}$ on $\partial B_{2}.$ We
let $\eta:(0,\infty)\to[0,1]$ be a $C^{1}$ function that satisfies
$\eta(t)\equiv 1$ for
$t\in\bigl{[}\nicefrac{{1}}{{\varepsilon^{2}}},\nicefrac{{1}}{{\varepsilon^{2}}}+\nicefrac{{b}}{{2}}\bigr{]},$
$\eta(t)\equiv 0$ when $t\geqslant\nicefrac{{1}}{{\varepsilon^{2}}}+b,$ and
$|\eta^{\prime}|\leqslant\frac{2}{b},$ and set

$\displaystyle\zeta(x):=-\Psi_{1}(x)\eta(|x-a_{2}^{\varepsilon}|)\,,$

where we recall that $a_{2}^{\varepsilon}$ is the center of
$B_{2}^{\varepsilon}.$ Then

$\displaystyle\nabla\zeta=-\eta(|x-a_{2}^{\varepsilon}|)\nabla\Psi_{1}(x)-\Psi_{1}(x)\eta^{\prime}(|x-a_{2}^{\varepsilon}|)\frac{x-a_{2}^{\varepsilon}}{|x-a_{2}^{\varepsilon}|},$

so that, pointwise, we have the bound

$\displaystyle|\nabla\zeta(x)|\leqslant\bigl{[}|\nabla\Psi_{1}(x)|+\frac{2}{b}|\Psi_{1}(x)|\bigr{]},$

with support in the set
$\nicefrac{{1}}{{\varepsilon^{2}}}\leqslant|x-a_{2}^{\varepsilon}|\leqslant\nicefrac{{1}}{{\varepsilon^{2}}}+b.$
Then, the energy of $\zeta$ is easily calculated:

$\displaystyle
F_{\varepsilon}(\zeta)\lesssim\int_{\nicefrac{{1}}{{\varepsilon^{2}}}\leqslant|x-a_{2}^{\varepsilon}|\leqslant\nicefrac{{1}}{{\varepsilon^{2}}}+b}\biggl{(}\frac{C}{b^{2}}|\Psi_{1}|^{2}+|\nabla\Psi_{1}|^{2}\biggr{)}\,dx$
(2.11)
$\displaystyle\quad=\biggl{|}\frac{C}{b^{2}}\int_{\nicefrac{{1}}{{\varepsilon^{2}}}\leqslant|x-a_{2}^{\varepsilon}|\leqslant\nicefrac{{1}}{{\varepsilon^{2}}}+b}|\Psi_{1}|^{2}\,dx+\int_{|x-a_{2}^{\varepsilon}|=\nicefrac{{1}}{{\varepsilon^{2}}}}\Psi_{1}(x)\frac{\partial\Psi_{1}(x)}{\partial\nu}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad\quad\quad-\int_{|x-a_{2}^{\varepsilon}|=\nicefrac{{1}}{{\varepsilon^{2}}}+b}\Psi_{1}(x)\frac{\partial\Psi_{1}(x)}{\partial\nu}\,d\mathcal{H}^{1}\biggr{|},$

where we plugged in the PDE satisfied by $\Psi_{1}$ and integrated by parts as
before; the signs in front of the boundary integrals reflect our choice that
the corresponding unit normals point towards $a_{2}^{\varepsilon}.$ Each of
these integrals can be explicitly estimated (see the proof of Lemma 2.4 for
details), and the triangle inequality then implies that

$|F_{\varepsilon}(\zeta)|\leqslant C\frac{K_{0}^{2}(2b)}{\varepsilon}\,.$

This completes the argument. ∎

In the next Lemma, we use Lemma 2.2 to control the boundary integral on the
right-hand side of (2.6).

######  Lemma 2.3.

For all $\varepsilon>0$ sufficiently small, we have the estimate

$\displaystyle\biggl{|}\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\biggl{(}\frac{\partial\Psi_{1}}{\partial\nu_{2}}+\frac{\partial
Z_{1}}{\partial\nu_{2}}\biggr{)}\biggr{|}\leqslant
C\frac{K_{0}^{2}(2b)}{\varepsilon}.$ (2.12)

###### Proof.

_Step 1._ By the triangle inequality,

$\displaystyle\biggl{|}\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\biggl{(}\frac{\partial\Psi_{1}}{\partial\nu_{2}}+\frac{\partial
Z_{1}}{\partial\nu_{2}}\biggr{)}\biggr{|}\leqslant\biggl{|}\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}\biggr{|}+\biggl{|}\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial
Z_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}\biggr{|}=:R_{1}+R_{2}.$

The previous lemma shows that the term $R_{1}$ is controlled by
$C\frac{K_{0}^{2}(2b)}{\varepsilon},$ so that the proof of the Lemma is
completed if we show the same bound for the term $R_{2}$.

_Step 2._ First we make the observation that the prescribed boundary
conditions on $Z_{1}$ imply that

$\displaystyle\int_{\partial B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial
Z_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}=\int_{\partial
B_{2}^{\varepsilon}}-Z_{1}\frac{\partial
Z_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}=-\int_{\partial\Omega_{\varepsilon}}Z_{1}\frac{\partial
Z_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}$

Then as $\Delta Z_{1}=Z_{1}$, integrating by parts we see that

$\displaystyle\int_{\partial\Omega^{\varepsilon}}Z_{1}\frac{\partial
Z_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}=\int_{\Omega^{\varepsilon}}|\nabla
Z_{1}|^{2}+|Z_{1}|^{2}\,dx=\|Z_{1}\|_{H^{1}(\Omega^{\varepsilon})}^{2}.$

Thus $R_{2}=\|Z_{1}\|_{H^{1}(\Omega^{\varepsilon})}^{2}$, and following the
estimates of Lemma 2.2 we thus have that

$\displaystyle R_{2}=\|Z_{1}\|_{H^{1}(\Omega^{\varepsilon})}^{2}\leqslant
C\frac{K_{0}(2b)^{2}}{\varepsilon}\,.$

∎

######  Lemma 2.4.

The off-diagonal terms from (2.7) have the asymptotic expansion

$\displaystyle\Biggl{|}\int_{\partial
B_{2}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{2}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\int_{\partial
B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}-\frac{2e^{-2b}\sqrt{\pi}}{\varepsilon}\Biggr{|}\leqslant
C\frac{K_{0}^{2}(2b)}{\varepsilon}\,.$ (2.13)

###### Proof.

_Step 1._ We proceed by a similar argument to Step 2 of Lemma 2.3. First we
note that

$\displaystyle\int_{\partial B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}=\int_{\partial
B_{1}^{\varepsilon}}-Z_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}=\int_{\partial\Omega_{\varepsilon}}-Z_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}$

Thus we may estimate this via

$\displaystyle\mathopen{}\mathclose{{}\left|\int_{\partial
B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}}\right|\leqslant\|Z_{2}\|_{H^{\nicefrac{{1}}{{2}}}(\partial\Omega)_{\varepsilon}}\Bigl{\|}\frac{\partial
Z_{1}}{\partial\nu}\Bigr{\|}_{H^{-\nicefrac{{1}}{{2}}}(\partial\Omega_{\varepsilon})}.$

We now utilise the fact that $\Delta Z_{i}=Z_{i}$ and the definition of the
$H^{\nicefrac{{1}}{{2}}}$ norm to conclude that
$\|Z_{2}\|_{H^{\nicefrac{{1}}{{2}}}(\partial\Omega_{\varepsilon})}=\|Z_{2}\|_{H^{1}(\Omega_{\varepsilon})}$,
and we use Proposition A.4 to conclude that $\Bigl{\|}\frac{\partial
Z_{1}}{\partial\nu}\Bigr{\|}_{H^{-\nicefrac{{1}}{{2}}}(\partial\Omega)}=\|Z_{1}\|_{H^{1}(\Omega_{\varepsilon})}$.
Thus by taking the estimations of the $H^{1}$ norms of $Z_{1},Z_{2}$ from
Lemma 2.2, we have that

$\displaystyle\mathopen{}\mathclose{{}\left|\int_{\partial
B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}}\right|\leqslant
C\frac{K_{0}(2b)^{2}}{\varepsilon}.$

_Step 2._ Towards evaluating the first two terms, first we write $\partial
B_{2}^{\varepsilon}=\partial B_{2}^{\varepsilon+}\cup\partial
B_{2}^{\varepsilon-},$ where the $\pm$ respectively denote the upper/lower
hemispheres (i.e., $y>-b-\frac{1}{\varepsilon^{2}}$ and
$y<-b-\frac{1}{\varepsilon^{2}}$ respectively). It is clear that the
contribution of $\partial B_{2}^{\varepsilon-}$ is exponentially small by
prior arguments, so we focus on the contribution of $\partial
B_{2}^{\varepsilon+}$ from the first two terms. We parameterize $\partial
B_{2}^{\varepsilon+}$ as a graph over the $x-$ axis:

$\displaystyle
y=-b-\frac{1}{\varepsilon^{2}}+\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}},\quad\quad|x|\leqslant\frac{1}{\varepsilon^{2}},$

where we choose the positive square root since we want the upper semicircle of
$\partial B_{2}^{\varepsilon}.$ We note that for $(x,y)\in\partial
B_{2}^{\varepsilon},$

$\displaystyle\nu_{2}=-\varepsilon^{2}\biggl{(}x,y+b+\frac{1}{\varepsilon^{2}}\biggr{)}.$

As a sign check, we note that at $x=0,y=-b$ the normal
$\nu_{2}=-\varepsilon^{2}(0,\frac{1}{\varepsilon^{2}})=(0,-1).$ Since

$\nabla\Psi_{1}(x,y)=\frac{1}{K_{0}\bigl{(}\frac{1}{\varepsilon^{2}}\bigr{)}}K_{0}^{\prime}\biggl{(}\sqrt{x^{2}+\bigl{(}y-b-\frac{1}{\varepsilon^{2}})^{2}\bigr{)}}\biggr{)}\frac{(x,y-b-\frac{1}{\varepsilon^{2}})}{\sqrt{x^{2}+\bigl{(}y-b-\frac{1}{\varepsilon^{2}}\bigr{)}^{2}}}\,,$

we arrive at

$\displaystyle\nu_{2}\cdot\nabla\Psi_{1}$
$\displaystyle=-\frac{\varepsilon^{2}}{K_{0}(\tfrac{1}{\varepsilon^{2}})}K_{0}^{\prime}\biggl{(}\biggl{\\{}x^{2}+\biggl{\\{}2b+\frac{2}{\varepsilon^{2}}-\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\biggr{\\}}^{2}\biggr{\\}}^{\nicefrac{{1}}{{2}}}\biggr{)}$
(2.14)
$\displaystyle\quad\quad\times\mathopen{}\mathclose{{}\left[\frac{x^{2}+\biggl{(}b+\frac{1}{\varepsilon^{2}}-\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\biggr{)}^{2}-(b+\frac{1}{\varepsilon^{2}})^{2}}{\biggl{\\{}x^{2}+\biggl{\\{}2b+\frac{2}{\varepsilon^{2}}-\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\biggr{\\}}^{2}\biggr{\\}}^{\nicefrac{{1}}{{2}}}}}\right],$

holding for all $(x,y)\in\partial B_{2}^{\varepsilon}.$ Note, with the
parametrization
$x\mapsto\bigl{(}x,-b-\frac{1}{\varepsilon^{2}}+\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\bigr{)},$
the speed of the curve $\partial B_{2}^{\varepsilon+}$ is
$\sqrt{1+\bigl{(}\frac{x}{\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}}\bigr{)}^{2}}=\bigl{(}\frac{1}{1-\varepsilon^{4}x^{2}}\bigr{)}^{\nicefrac{{1}}{{2}}}.$

_Step 3._ We will split the integral as $[0,\frac{M}{\varepsilon}],$ and
between $[\frac{M}{\varepsilon},\frac{1}{\varepsilon^{2}}],$ for some
$M\leqslant\frac{\sqrt{3}}{2}\frac{1}{\varepsilon}<\frac{1}{\varepsilon},$ to
be fixed later. We compute each of these contributions separately. For the
first integral in (2.13) we obtain

$\int_{\partial
B_{2}^{\varepsilon+}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}=\biggl{(}2\int_{0}^{\frac{M}{\varepsilon}}+\int_{|x|\geqslant\frac{M}{\varepsilon}}\biggr{)}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}\,.$
(2.15)

Towards computing it we note by the binomial theorem that for
$|x|\leqslant\frac{M}{\varepsilon}$,

$\displaystyle\biggl{\\{}x^{2}+\biggl{\\{}2b+\frac{2}{\varepsilon^{2}}-\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\biggr{\\}}^{2}\biggr{\\}}^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\biggl{\\{}\frac{1}{\varepsilon^{4}}+\frac{4}{\varepsilon^{4}}(1+b\varepsilon^{2})\Bigl{(}1+b\varepsilon^{2}-\sqrt{1-\varepsilon^{4}x^{2}}\Bigr{)}\biggr{\\}}^{\nicefrac{{1}}{{2}}}$
$\displaystyle=:\biggl{\\{}\frac{1}{\varepsilon^{4}}+C_{M}(x)\biggr{\\}}^{\nicefrac{{1}}{{2}}}$

where the function $C_{M}$ satisfies

  1. (i)

$C_{M}(0)=\frac{4b}{\varepsilon^{2}}(1+b\varepsilon^{2}),$

  2. (ii)

$C_{M}(\frac{\sqrt{3}}{2\varepsilon})\geqslant\frac{3}{\varepsilon^{4}},$ so
that $\sqrt{\varepsilon^{-4}+C_{M}(x)}\geqslant\frac{2}{\varepsilon^{2}}$ for
$x=\frac{\sqrt{3}}{2\varepsilon},$ and

  3. (iii)

$C_{M}(x)$ is an increasing function for all
$x\in[0,\frac{1}{\varepsilon^{2}}],$ so that if
$x\geqslant\frac{\sqrt{3}}{2\varepsilon},$ then the preceding quantity, i.e.,
$\sqrt{\varepsilon^{-4}+C_{M}(x)}\geqslant\frac{2}{\varepsilon^{2}}.$

In addition, we have

$\displaystyle\frac{x^{2}+\biggl{(}b+\frac{1}{\varepsilon^{2}}-\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\biggr{)}^{2}-(b+\frac{1}{\varepsilon^{2}})^{2}}{\biggl{\\{}x^{2}+\biggl{\\{}2b+\frac{2}{\varepsilon^{2}}-\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\biggr{\\}}^{2}\biggr{\\}}^{\nicefrac{{1}}{{2}}}}\frac{1}{\sqrt{1-\varepsilon^{4}x^{2}}}$
$\displaystyle\quad\quad=\frac{-\frac{1}{\varepsilon^{4}}+\frac{2}{\varepsilon^{4}}\bigl{(}1-(1+b\varepsilon^{2})\sqrt{1-\varepsilon^{4}x^{2}}\bigr{)}}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}\frac{1}{\sqrt{1-\varepsilon^{4}x^{2}}}\,.$

It follows then that the first term in (2.15) contributes (see (2.14))

$\displaystyle
2\int_{0}^{\frac{M}{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}$
$\displaystyle\quad=-2\varepsilon^{2}\int_{0}^{\frac{M}{\varepsilon}}\frac{1}{K_{0}(\frac{1}{\varepsilon^{2}})}K_{0}^{\prime}\biggl{(}\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}\biggr{)}\frac{-\frac{1}{\varepsilon^{4}}+\frac{2}{\varepsilon^{4}}\bigl{(}1-(1+b\varepsilon^{2})\sqrt{1-\varepsilon^{4}x^{2}}\bigr{)}}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}\frac{1}{\sqrt{1-\varepsilon^{4}x^{2}}}\,dx$
$\displaystyle\quad=-2\int_{0}^{\frac{M}{\varepsilon}}\frac{1}{K_{0}(\frac{1}{\varepsilon^{2}})}K_{0}^{\prime}\biggl{(}\frac{1}{\varepsilon^{2}}\sqrt{1+\varepsilon^{4}C_{M}(x)}\biggr{)}\frac{-1+2\bigl{(}1-(1+b\varepsilon^{2})\sqrt{1-\varepsilon^{4}x^{2}}\bigr{)}}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}\frac{1}{\sqrt{1-\varepsilon^{4}x^{2}}}\,dx$
$\displaystyle\quad=2\int_{0}^{\frac{M}{\varepsilon}}\frac{1}{K_{0}(\frac{1}{\varepsilon^{2}})}K_{0}^{\prime}\biggl{(}\frac{1}{\varepsilon^{2}}\sqrt{1+\varepsilon^{4}C_{M}(x)}\biggr{)}\,dx+R_{1},$

where the remainder $R_{1}$ satisfies

$|R_{1}|\leqslant
2\varepsilon^{2}\Biggl{|}\int_{0}^{\frac{M}{\varepsilon}}\frac{1}{K_{0}(\frac{1}{\varepsilon^{2}})}K_{0}^{\prime}\biggl{(}\frac{1}{\varepsilon^{2}}\sqrt{1+\varepsilon^{4}C_{M}(x)}\biggr{)}\,dx\Biggr{|}\,.$

Finally, by Lemma C.1, we obtain

$\displaystyle
2\int_{0}^{\frac{M}{\varepsilon}}\frac{1}{K_{0}(\frac{1}{\varepsilon^{2}})}K_{0}^{\prime}\biggl{(}\frac{1}{\varepsilon^{2}}\sqrt{1+\varepsilon^{4}C_{M}(x)}\biggr{)}\,dx$
$\displaystyle\quad\approx
2\int_{0}^{\frac{M}{\varepsilon}}\exp\Bigl{(}\frac{1}{\varepsilon^{2}}\bigl{(}1-\sqrt{1+\varepsilon^{4}C_{M}(x)}\bigr{)}\Bigr{)}\sqrt{\frac{1}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}}\,dx\,.$

Inserting the definition of
$C_{M}(x)=\frac{4}{\varepsilon^{4}}(1+b\varepsilon^{2})\bigl{(}1+b\varepsilon^{2}-\sqrt{1-\varepsilon^{4}x^{2}}\bigr{)},$
and making the change of variables $s=\varepsilon^{2}x,$ we find that the last
integral simplifies to

$\displaystyle\frac{2}{\varepsilon^{2}}\int_{0}^{M\varepsilon}\\!\\!\exp\Bigl{(}\frac{-\frac{4}{\varepsilon^{2}}(1+b\varepsilon^{2})(1+b\varepsilon^{2}-\sqrt{1-s^{2}})}{1+\sqrt{1+4(1+b\varepsilon^{2})(1+b\varepsilon^{2}-\sqrt{1-s^{2}})}}\Bigr{)}\Bigl{(}1{+}4(1{+}b\varepsilon^{2})(1{+}b\varepsilon^{2}{-}\sqrt{1{-}s^{2}})\Bigr{)}^{-\frac{1}{4}}\,ds$
$\displaystyle\quad\approx\frac{2}{\varepsilon^{2}}\int_{0}^{M\varepsilon}\exp\Bigl{(}-\frac{2}{\varepsilon^{2}}(b\varepsilon^{2}+\frac{1}{2}s^{2})\Bigr{)}\,ds$
$\displaystyle\quad=\frac{2e^{-2b}}{\varepsilon^{2}}\int_{0}^{M\varepsilon}\exp\Bigl{(}-\frac{s^{2}}{\varepsilon^{2}}\Bigr{)}\,ds$
$\displaystyle\quad=\frac{2e^{-2b}}{\varepsilon}\int_{0}^{M}\exp(-s^{2})\,ds=\frac{e^{-2b}\sqrt{\pi}}{\varepsilon}\Bigl{[}1-\frac{2}{\sqrt{\pi}}\int_{M}^{\infty}\exp(-s^{2})\,ds\Bigr{]}=\frac{e^{-2b}\sqrt{\pi}}{\varepsilon}\bigl{(}1-\mathrm{erfc}(M)\bigr{)}\,.$

In each of the preceding two displays, the $\approx$ (approximate) sign means
that the left-hand and right-hand sides differ by $O(\varepsilon).$ To
conclude the computation of the leading term in (2.15), we observe that the
complimentary error function $\mathrm{erfc}$ satisfies the asymptotics

$\Bigl{|}\mathrm{erfc}(M)-\frac{e^{-M^{2}}}{M\sqrt{\pi}}\Bigr{|}\leqslant
CM^{-3}e^{-M^{2}}\,.$

At this point, we must choose $M$ so that $CM^{-3}e^{-M^{2}}\ll 1$ as
$\varepsilon\to 0.$ This is, for example, guaranteed with the choice111the
precise prefactor is not important, but is chosen to simplify the arithmetic
in our computation of the tails $M=\frac{\sqrt{3}}{2}\frac{1}{\varepsilon},$
so that combining the preceding displays we find

$\Biggl{|}2\int_{0}^{\frac{M}{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}-\frac{\sqrt{\pi}e^{-2b}}{\varepsilon}\Biggr{|}\leqslant
C\varepsilon+\frac{Ce^{-2b}}{\varepsilon^{4}}e^{-\frac{C}{\varepsilon^{2}}}\leqslant
C\varepsilon\,.$

For the tail term in (2.15), noting from the properties of the function
$C_{M},$ that
$\sqrt{\varepsilon^{-4}+C_{M}(\frac{\sqrt{3}}{2\varepsilon})}\geqslant\frac{2}{\varepsilon^{2}},$
and that $C_{M}$ is increasing, we find

$\biggl{|}\int_{|x|\geqslant\frac{M}{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}\biggr{|}\leqslant
C\int_{\frac{M}{\varepsilon}}^{\frac{1}{\varepsilon^{2}}}\exp\bigl{(}\frac{1}{\varepsilon^{2}}(1-\sqrt{1+\varepsilon^{4}C_{M}(x)})\bigr{)}\,dx\leqslant\frac{1}{\varepsilon^{2}}\exp\biggl{(}-\frac{1}{\varepsilon^{2}}\biggr{)}\leqslant
C\varepsilon\,,$

for all $\varepsilon$ small enough.

_Step 4._ The last step in the proof is to evaluate the asymptotics of the
term

$\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{2}}{\partial\nu_{2}}\,d\mathcal{H}^{1}\,.$

This is easier than Step 3, since

$\frac{\partial\Psi_{2}}{\partial\nu_{2}}=-\frac{1}{K_{0}\bigl{(}\frac{1}{\varepsilon^{2}}\bigr{)}}K_{0}^{\prime}\Bigl{(}\frac{1}{\varepsilon^{2}}\Bigr{)}\,.$

Consequently, we find

$\displaystyle\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{2}}{\partial\nu_{2}}\,d\mathcal{H}^{1}$
$\displaystyle\quad=-2\frac{1}{K_{0}\bigl{(}\frac{1}{\varepsilon^{2}}\bigr{)}}K_{0}^{\prime}\Bigl{(}\frac{1}{\varepsilon^{2}}\Bigr{)}\int_{0}^{\frac{1}{\varepsilon^{2}}}\frac{1}{K_{0}\bigl{(}\frac{1}{\varepsilon^{2}}\bigr{)}}K_{0}\Bigl{(}\frac{1}{\varepsilon^{2}}\sqrt{1+\varepsilon^{4}C_{M}(x)}\Bigr{)}\frac{1}{\sqrt{1-\varepsilon^{4}x^{2}}}\,dx+\frac{1}{\varepsilon^{2}}\exp\Bigl{(}-\frac{1}{\varepsilon^{2}}\Bigr{)}$
$\displaystyle\approx\frac{e^{-2b}\sqrt{\pi}}{\varepsilon}+C\varepsilon+\frac{1}{\varepsilon^{2}}\exp\Bigl{(}-\frac{1}{\varepsilon^{2}}\Bigr{)}\,.$

The proof of the proposition is completed by combining Steps 1 through 4. ∎

###### Proof of Theorem 1.1.

Theorem 1.1 easily follows by combining (2.5), and Lemmas 2.1, 2.3 and 2.4. ∎

###  Interaction energies for O(1) separation between particles

In this section we use formal asymptotics to compute the energy of interaction
between two particles when $b=O(1)$ and $\varepsilon\to 0^{+}$. This amounts
to computing various terms in (2.6) and (2.7). Since the problem is
rotationally invariant, in this section we find it convenient to orient the
particles horizontally (Fig. 2), rather than vertically. Given
$b,\varepsilon>0,$ consider two disk-like particles $B_{1}^{\varepsilon}$ and
$B_{2}^{\varepsilon}$ of radius $1/\varepsilon^{2}$ where the first particle
is centered at the origin and the distance between the particles is equal $2b$
as shown in Fig. 2.

Figure 2: Geometry of the problem

We begin by introducing polar coordinates associated with the center of the
particle $B_{1}^{\varepsilon}$ so that

$x=-\varepsilon^{-2}-b+\rho\cos{\theta},\quad y=\rho\sin{\theta},$

then

$\partial
B_{1}^{\varepsilon}=\mathopen{}\mathclose{{}\left\\{(\rho,\theta)\mathopen{}\mathclose{{}\left|\,\rho=\varepsilon^{-2},\,\theta\in[0,2\pi)}\right.}\right\\}$

and

$\partial
B_{2}^{\varepsilon}=\mathopen{}\mathclose{{}\left\\{(\rho,\theta)\mathopen{}\mathclose{{}\left|\,\rho^{2}-4(\varepsilon^{-2}+b)\rho\cos{\theta}+4{(\varepsilon^{-2}+b)}^{2}-\varepsilon^{-4}=0,\,\theta\in[0,2\pi)}\right.}\right\\}.$

Now suppose that $\varepsilon\ll 1$ and let $t:=\varepsilon^{-1}\,\theta.$
Then, if $t=O(1),$ we have
$\theta=O\mathopen{}\mathclose{{}\left(\varepsilon}\right)$ so that

$\cos{\theta}=1-\frac{\varepsilon^{2}t^{2}}{2}+O\mathopen{}\mathclose{{}\left(\varepsilon^{4}}\right)$

and the equation for $\partial B_{2}^{\varepsilon}$ is

$\rho^{2}-4(\varepsilon^{-2}+b)\rho\mathopen{}\mathclose{{}\left(1-\frac{\varepsilon^{2}t^{2}}{2}}\right)+4{(\varepsilon^{-2}+b)}^{2}-\varepsilon^{-4}=0,$

up to the order $O\mathopen{}\mathclose{{}\left(\varepsilon}\right).$ Solving
this equation for $\rho$, gives an asymptotic expression for the boundary of
the right disk, i.e.,

$\rho|_{\partial B_{2}^{\varepsilon}}=\varepsilon^{-2}+2b+t^{2},$ (2.16)

valid up to $O\mathopen{}\mathclose{{}\left(\varepsilon}\right),$ while the
boundary of the left disk is given by

$\rho|_{\partial B_{1}^{\varepsilon}}=\varepsilon^{-2}.$

We now solve the problem (2.2) for the left disk

$\mathopen{}\mathclose{{}\left\\{\begin{aligned}
&-\Delta\Psi_{1}+\Psi_{1}=0\mbox{ in }\mathbb{R}^{2}\backslash
B_{1}^{\varepsilon},\\\ &\Psi_{1}|_{\partial
B_{1}^{\varepsilon}}=1.\end{aligned}}\right.$ (2.17)

Using the radially symmetric ansatz $\Psi_{1}=\Psi_{1}(\rho)$ and the fact
that $\rho>\varepsilon^{-2}\gg 1,$ we find that

$\Psi_{1}(\rho)\sim e^{\varepsilon^{-2}-\rho}.$ (2.18)

to leading order. Note that this result matches the expression (2.1) for small
$\varepsilon.$ Therefore

$\Psi_{1}|_{\partial
B_{2}^{\varepsilon}}=\Psi_{1}\mathopen{}\mathclose{{}\left(\varepsilon^{-2}+2b+t^{2}}\right)=e^{-2b-t^{2}}$
(2.19)

and this expression decays exponentially fast in $t$ as one moves away from
the point on $\partial B_{2}^{\varepsilon}$ that is the closest to $\partial
B_{1}^{\varepsilon}.$

Now we can solve the problem (2.3). We have

$\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}-Z_{1,\rho\rho}-\frac{1}{\rho}Z_{1,\rho}-\frac{1}{\varepsilon^{2}\rho^{2}}Z_{1,tt}+Z_{1}=0,&(\rho,t)\in\Omega_{\varepsilon}\\\
Z_{1}\mathopen{}\mathclose{{}\left(\rho,\varepsilon
t}\right)=0,&\rho=\varepsilon^{-2},\\\
Z_{1}\mathopen{}\mathclose{{}\left(\rho,\varepsilon
t}\right)=-e^{-2b-t^{2}},&(\rho,t)\in\partial
B_{2}^{\varepsilon}.\end{array}}\right.$ (2.20)

Because $\rho>\varepsilon^{-2},$ to leading order the equation in (2.20) takes
the form

$-Z_{1,\rho\rho}+Z_{1}=0,$

with the general solution

$Z_{1}=C_{1}e^{-\rho}+C_{2}e^{\rho}.$

Substituting this $Z_{1}$ into the boundary conditions from (2.20) gives

$\mathopen{}\mathclose{{}\left\\{\begin{array}[]{l}C_{1}e^{-\varepsilon^{-2}}+C_{2}e^{\varepsilon^{-2}}=0,\\\
C_{1}e^{-\varepsilon^{-2}-2b-t^{2}}+C_{2}e^{\varepsilon^{-2}+2b+t^{2}}=-e^{-2b-t^{2}},\end{array}}\right.$

so that

$C_{1}=\frac{e^{\varepsilon^{-2}}}{e^{4b+2t^{2}}-1},\quad
C_{2}=-\frac{e^{-\varepsilon^{-2}}}{e^{4b+2t^{2}}-1},$

hence

$Z_{1}=\frac{2\sinh{(\varepsilon^{-2}-\rho)}}{e^{4b+2t^{2}}-1}$ (2.21)

solves (2.20) to leading order in $\varepsilon.$ It then follows that

$\mathopen{}\mathclose{{}\left.\frac{\partial
Z_{1}}{\partial\nu_{1}}}\right|_{\partial
B_{1}^{\varepsilon}}=-\mathopen{}\mathclose{{}\left.\frac{\partial
Z_{1}}{\partial\rho}}\right|_{\partial
B_{1}^{\varepsilon}}=\frac{2}{e^{4b+2t^{2}}-1}.$ (2.22)

Next, we observe that

$\Psi_{1}\mathopen{}\mathclose{{}\left(\rho(\theta)}\right)=\Psi_{2}\mathopen{}\mathclose{{}\left(\bar{\rho}(\bar{\theta})}\right)$

(cf. Fig. 3).

Figure 3: Setup for the calculation of $z$

Recalling that $\theta=\varepsilon t,$ we set
$\bar{\theta}=\varepsilon\bar{t},$ use the law of sines and (2.16) to find
that

$\frac{\varepsilon^{-2}+2b+t^{2}}{\sin{(\varepsilon\bar{t})}}=\frac{1}{\varepsilon^{2}\sin{(\varepsilon
t)}},$

hence

$\mathopen{}\mathclose{{}\left(\varepsilon^{-2}+2b+t^{2}}\right)t=\varepsilon^{-2}\bar{t},$

because $\varepsilon^{-2}t,\varepsilon^{-2}\bar{t}\ll 1.$ We conclude that

$\bar{t}=t+O\mathopen{}\mathclose{{}\left(\varepsilon^{2}}\right)$

and

$\Psi_{2}|_{\partial
B_{1}^{\varepsilon}}=\Psi_{1}\mathopen{}\mathclose{{}\left(\rho\mathopen{}\mathclose{{}\left(\varepsilon
t}\right)}\right)+O\mathopen{}\mathclose{{}\left(\varepsilon^{2}}\right)=e^{-2b-t^{2}}+O\mathopen{}\mathclose{{}\left(\varepsilon^{2}}\right).$
(2.23)

We are now ready to determine the $\varepsilon\to 0^{+}$ asymptotics of terms
that appear in (2.6-2.7) in a series of six Steps:

Step (1). Using (2.18), we obtain

$\int_{\partial
B_{1}}\frac{\partial\Psi_{1}}{\partial\nu_{1}}\sim{\frac{2\pi}{\varepsilon^{2}}}$

Step(2). From (2.22) and (2.23), we deduce

$\int_{\partial B_{1}^{\varepsilon}}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}=\varepsilon^{-2}\int_{-\pi}^{\pi}\Psi_{2}\frac{\partial
Z_{1}}{\partial\nu_{1}}\,d\theta=\varepsilon^{-1}\int_{-\pi/\varepsilon}^{\pi/\varepsilon}\Psi_{2}\frac{\partial{Z}_{1}}{\partial\nu_{1}}\,dt\\\
\sim\frac{2}{\varepsilon}\int_{-\infty}^{\infty}\frac{e^{-6b}e^{-3t^{2}}}{1-e^{-4b}e^{-2t^{2}}}\,dt\sim\frac{\sqrt{2\pi}}{\varepsilon}\Theta_{3}\mathopen{}\mathclose{{}\left(e^{-2b}}\right),$
(2.24)

where

$\Theta_{k}(x):=\sqrt{\frac{2}{\pi}}\int_{-\infty}^{\infty}\frac{x^{k}e^{-kt^{2}}}{1-x^{2}e^{-2t^{2}}}\,dt,$
(2.25)

where $x\in(0,1)$ and $k\in\mathbb{N}$. Note that
${\Theta_{2}(x)=\mathrm{Li}_{\frac{1}{2}}\,\mathopen{}\mathclose{{}\left(x^{2}}\right)},$
where $\mathrm{Li}_{\frac{1}{2}}\,$ is the polylogarithmic function and

$\lim_{x\to 0}{\frac{\Theta_{k}(x)}{x^{k}}}=\sqrt{\frac{2}{k}}.$ (2.26)

Step(3). Using (2.19), we have

$\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{2}}{\partial\nu_{2}}\,d\mathcal{H}^{1}=\varepsilon^{-2}\frac{\partial\Psi_{2}}{\partial\nu_{2}}(\varepsilon^{-2})\int_{-\pi}^{\pi}\Psi_{1}\,d\theta\\\
=\frac{1}{\varepsilon^{2}}\frac{\partial\Psi_{2}}{\partial\nu_{2}}(\varepsilon^{-2})\int_{-\pi/\varepsilon}^{\pi/\varepsilon}\Psi_{1}\,dt\sim\frac{e^{-2b}}{\varepsilon^{2}}\int_{-\infty}^{\infty}e^{-t^{2}}\,dt=\frac{\sqrt{\pi}e^{-2b}}{\varepsilon}.$
(2.27)

Figure 4: Setup for the calculation in (5).

Step (4). Using Fig. 4 and the law of cosines, we have

$\rho^{2}=4{(\varepsilon^{-2}+b)}^{2}+\rho_{1}^{2}-4\rho_{1}(\varepsilon^{-2}+b)\cos{\theta_{1}},$
(2.28)

then $\Psi_{1}$ can be written as

$\bar{\Psi}_{1}(\rho_{1},\theta_{1})=\Psi_{1}(\rho(\rho_{1},\theta_{1}),\theta(\rho_{1},\theta_{1}))\sim
e^{\varepsilon^{-2}-\sqrt{4{(\varepsilon^{-2}+b)}^{2}+\rho_{1}^{2}-4\rho_{1}(\varepsilon^{-2}+b)\cos{\theta_{1}}}}$

so that

$\frac{\partial\Psi_{1}}{\partial\nu_{2}}=-\mathopen{}\mathclose{{}\left.\frac{\partial\bar{\Psi}_{1}}{\partial\rho_{1}}}\right|_{\rho_{1}=\varepsilon^{-2}}\\\
\sim\frac{2(1+b\varepsilon^{2})\cos{\theta_{1}}-1}{\sqrt{4{(1+b\varepsilon^{2})}^{2}+1-4(1+b\varepsilon^{2})\cos{\theta_{1}}}}e^{\varepsilon^{-2}(1-\sqrt{4{(1+b\varepsilon^{2})}^{2}+1-4(1+b\varepsilon^{2})\cos{\theta_{1}}})}$

and

$\frac{\partial\Psi_{1}}{\partial\nu_{2}}\sim-e^{-2b-t_{1}^{2}}.$

Here the last step follows because

$1-2\cos{\theta_{1}}\sim-1+\varepsilon^{2}t_{1}^{2}$

and

$4{(\varepsilon^{-2}+b)}^{2}+\varepsilon^{-4}-4\varepsilon^{-2}(\varepsilon^{-2}+b)\cos{\theta_{1}}\sim\varepsilon^{-2}+2b+t_{1}^{2},$

where $t_{1}:=\theta_{1}/\varepsilon.$ Therefore

$\int_{\partial
B_{2}^{\varepsilon}}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}\sim-\frac{e^{-2b}}{\varepsilon}\int_{-\pi/\varepsilon}^{\pi/\varepsilon}e^{-t_{1}^{2}}\,dt_{1}\sim-\frac{\sqrt{\pi}}{\varepsilon}e^{-2b}$
(2.29)

Step(5). Using (2.19), we find

$\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial\Psi_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}\sim-\frac{e^{-4b}}{\varepsilon}\int_{-\pi/\varepsilon}^{\pi/\varepsilon}e^{-2t^{2}}\,dt\sim-\frac{\sqrt{\pi/2}}{\varepsilon}e^{-4b}$
(2.30)

Step(6). The last term that we need to estimate is $\int_{\partial
B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial Z_{1}}{\partial\nu_{2}}.$ From
(2.21), we have that

$\bar{Z}_{1}(\rho_{1},t_{1})=Z_{1}\mathopen{}\mathclose{{}\left(\rho(\rho_{1},t_{1}),t(\rho_{1},t_{1})}\right)\sim\frac{2\sinh{(\varepsilon^{-2}-\rho(\rho_{1},t_{1}))}}{e^{4b+2t(\rho_{1},t_{1})^{2}}-1},$

where $t_{1}=\theta_{1}/\varepsilon$ and $(\rho_{1},\theta_{1})$ are as in
Fig. 4. Further,

$\rho(\rho_{1},t_{1})=\sqrt{4{(\varepsilon^{-2}+b)}^{2}+\rho_{1}^{2}-4\rho_{1}(\varepsilon^{-2}+b)\cos{\theta_{1}}}$

per (2.28) and

$2(\varepsilon^{-2}+b)\sin{\theta}=\rho_{1}\sin{(\theta+\theta_{1}),}$

by the law of sines, hence

$\frac{\partial
Z_{1}}{\partial\nu_{2}}=-\mathopen{}\mathclose{{}\left.\frac{\partial\bar{Z}_{1}}{\partial\rho_{1}}}\right|_{\rho_{1}=\varepsilon^{-2}}\sim-\frac{2e^{-4b-2t_{1}^{2}}}{1-e^{-4b-2t_{1}^{2}}}\cosh{(2b+t_{1}^{2})}.$

It follows that

$\displaystyle\int_{\partial B_{2}^{\varepsilon}}\Psi_{1}\frac{\partial
Z_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}$
$\displaystyle\sim-\frac{2}{\varepsilon}\int_{-\pi/\varepsilon}^{\pi/\varepsilon}\frac{e^{-6b-3t_{1}^{2}}}{1-e^{-4b-2t_{1}^{2}}}\cosh{(2b+t_{1}^{2})}$
(2.31)
$\displaystyle\sim-\frac{\sqrt{\pi/2}}{\varepsilon}\mathopen{}\mathclose{{}\left(\mathrm{Li}_{\frac{1}{2}}\,\mathopen{}\mathclose{{}\left(e^{-4b}}\right)+\Theta_{4}\mathopen{}\mathclose{{}\left(e^{-2b}}\right)}\right).$
(2.32)

Now collecting the energy contributions in Steps (1) through (6) and using
(2.6-2.7), we find that

$\kappa_{\varepsilon}=g_{1}^{2}{\|\Psi_{1}+z_{1}\|}^{2}+g_{2}^{2}{\|\Psi_{2}+z_{2}\|}^{2}+2g_{1}g_{2}\langle\Psi_{1}+z_{1},\Psi_{2}+z_{2}\rangle\\\
=\frac{2\pi}{\varepsilon^{2}}\mathopen{}\mathclose{{}\left(g_{1}^{2}+g_{2}^{2}}\right)+\frac{\sqrt{\pi/2}}{\varepsilon}\mathopen{}\mathclose{{}\left[\mathrm{Li}_{\frac{1}{2}}\,\mathopen{}\mathclose{{}\left(e^{-4b}}\right)+\Theta_{4}\mathopen{}\mathclose{{}\left(e^{-2b}}\right)+e^{-4b}}\right](g_{1}^{2}+g_{2}^{2})\\\
-\frac{4\sqrt{\pi}}{\varepsilon}\mathopen{}\mathclose{{}\left[e^{-2b}+\frac{1}{\sqrt{2}}\Theta_{3}\mathopen{}\mathclose{{}\left(e^{-2b}}\right)}\right]g_{1}g_{2}+O(1)\,.$
(2.33)

Note that, using (2.26), this expression reduces to

$\kappa_{\varepsilon}=\frac{2\pi}{\varepsilon^{2}}\mathopen{}\mathclose{{}\left(g_{1}^{2}+g_{2}^{2}}\right)+\frac{\sqrt{\pi/2}}{\varepsilon}\mathopen{}\mathclose{{}\left[2e^{-4b}+\frac{1}{\sqrt{2}}e^{-8b}}\right](g_{1}^{2}+g_{2}^{2})\\\
-\frac{4\sqrt{\pi}}{\varepsilon}\mathopen{}\mathclose{{}\left[e^{-2b}+\frac{1}{\sqrt{3}}e^{-6b}}\right]g_{1}g_{2}+O(1),$

when $b\gg 1,$ matching the asymptotics of $\kappa_{\varepsilon}$ established
in Theorem 1.1.

######  Remark 2.5.

Consider a single particle $B_{1/\varepsilon^{2}}$ of the radius
$1/\varepsilon^{2}$ centered at the origin and let $k(T)=2$ in (1.1) so that

$W(u)=2{|u|}^{2}-2{|u|}^{4}+{|u|}^{6}.$

We can use asymptotic arguments developed in this section to compare the rates
of decay of solutions to the nonlinear

$\displaystyle\Delta u$ $\displaystyle=\frac{1}{4}\nabla_{u}W(u),\quad u\in
H^{1}\mathopen{}\mathclose{{}\left(\mathbb{R}^{2}\setminus\overline{B}_{1/\varepsilon^{2}}}\right)$
(2.34) $\displaystyle u$ $\displaystyle=1\qquad\mbox{ on }\partial
B_{1/\varepsilon^{2}}\,.$

and linear

$\displaystyle\Delta v$ $\displaystyle=v,\quad v\in
H^{1}\mathopen{}\mathclose{{}\left(\mathbb{R}^{2}\setminus\overline{B}_{1/\varepsilon^{2}}}\right)$
(2.35) $\displaystyle v$ $\displaystyle=1\qquad\mbox{ on }\partial
B_{1/\varepsilon^{2}}$

problems when $\varepsilon$ is small. Indeed, suppose that $(\rho,\theta)$ are
polar coordinates associated with the center of a particle
$B_{1/\varepsilon^{2}}$ and $\varepsilon\ll 1$. From the proof of Lemma B.1
for $m=0$ and (2.18), we have that the solution of (2.35) is

$v=\frac{K_{0}(\rho)}{K_{0}\big{(}\frac{1}{\varepsilon^{2}}\big{)}}\sim
e^{1/\varepsilon^{2}-\rho},$ (2.36)

because $\rho>\frac{1}{\varepsilon^{2}}$ and $\varepsilon\ll 1.$ Assuming that
$u=u(\rho)$ in the same asymptotic regime $u$ satisfies

$\displaystyle u_{\rho\rho}$ $\displaystyle=\frac{1}{4}\nabla_{u}W(u),$ (2.37)
$\displaystyle
u\mathopen{}\mathclose{{}\left(\frac{1}{\varepsilon^{2}}}\right)$
$\displaystyle=1,\quad\lim_{\rho\to\infty}u=0\,.$

to leading order in $\varepsilon.$ This problem has an explicit solution

$u=\frac{\sqrt{2}}{\sqrt{1+\cosh{2\mathopen{}\mathclose{{}\left(\rho-\frac{1}{\varepsilon^{2}}}\right)}+\sqrt{2}\sinh{2\mathopen{}\mathclose{{}\left(\rho-\frac{1}{\varepsilon^{2}}}\right)}}}\sim\frac{2}{\sqrt{1+\sqrt{2}}}e^{1/\varepsilon^{2}-\rho},$

when $\rho\gg\frac{1}{\varepsilon^{2}}$ and it would be reasonable to expect
that the interaction between two well-separated particles would be the same to
leading order, up to a multiplicative constant.

It is possible to quantify the tail behaviors of $u$ and $v$ rigorously for
general boundary conditions. At its heart is a convexity argument. We begin by
noticing that with $W$ as in (1.1) and $k(T)>\frac{4}{3},$ the function $W$ is
uniformly convex and $C^{2},$ so that for any $p,q\in\mathbb{R}^{2},$ we have

$\bigl{(}\nabla W(p)-\nabla W(q)\bigr{)}\cdot(p-q)\geqslant c_{0}|p-q|^{2}\,,$
(2.38)

for some $c_{0}>0$, that only depends on $k(T)>\frac{4}{3}.$ Then we have,

######  Proposition 2.6.

Let $v_{\varepsilon}$ denote the unique
$H^{1}(\mathbb{R}^{2}\setminus\overline{\widehat{B}^{\varepsilon}_{1}\cup\widehat{B}^{\varepsilon}_{2}})$
solution to the nonlinear system (1.3), and let $u_{\varepsilon}$ the
corresponding unique solution to the linear system (1.4), both with the same
boundary condition. Then, we claim that there exists $\alpha,C_{0}>0$ that
only depend on $c_{0},$ from (2.38) such that

$\int_{|x|\geqslant C_{0}|\log\varepsilon|}e^{2\alpha
H_{\varepsilon}(x)}|u_{\varepsilon}-v_{\varepsilon}|^{2}\,dx\lesssim
C\varepsilon^{4}\,,$ (2.39)

where

$H_{\varepsilon}(x)=\frac{1}{\varepsilon^{2}}\mathrm{dist}(x,\overline{\widehat{B}^{\varepsilon}_{1}\cup\widehat{B}^{\varepsilon}_{2}}).$
(2.40)

###### Proof of the Proposition 2.6.

Let $\eta$ be a positive smooth compactly supported test function that will be
subsequently chosen. We compute that

$\Delta
u_{\varepsilon}-\frac{1}{2\varepsilon^{4}}\nabla_{u}W(u_{\varepsilon})=\Delta(u_{\varepsilon}-v_{\varepsilon})+\frac{1}{2\varepsilon^{4}}\bigl{(}\nabla_{u}W(v_{\varepsilon})-\nabla_{u}W(u_{\varepsilon})\bigr{)}\,,$

so that taking the dot product of both sides with
$\eta^{2}(v_{\varepsilon}-u_{\varepsilon}),$ integrating on
$\mathbb{R}^{2}\setminus\overline{\widehat{B}^{\varepsilon}_{1}\cup\widehat{B}^{\varepsilon}_{2}}$,
using the uniform convexity of $W,$ integrating by parts, and using that $u=v$
on the boundary, yields

$\displaystyle
C\int|\nabla\eta|^{2}|u-v|^{2}+\int\eta^{2}(v_{\varepsilon}-u_{\varepsilon})\cdot\biggl{(}\Delta
u_{\varepsilon}-\frac{1}{2\varepsilon^{4}}\nabla_{u}W(u_{\varepsilon})\biggr{)}\,dx$
$\displaystyle\qquad\qquad\geqslant\frac{1}{2}\int\eta^{2}|\nabla
u_{\varepsilon}-\nabla
v_{\varepsilon}|^{2}+\frac{c_{0}}{2\varepsilon^{4}}\int\eta^{2}|v_{\varepsilon}-u_{\varepsilon}|^{2}\,dx\,.$

As $|u_{\varepsilon}|\leqslant 1$, we estimate

$\biggl{|}\Delta
u_{\varepsilon}-\frac{1}{2\varepsilon^{4}}\nabla_{u}W(u_{\varepsilon})\biggr{|}=\frac{1}{2\varepsilon^{4}}\bigl{|}\bigl{(}8|u_{\varepsilon}|^{2}-6|u_{\varepsilon}|^{4}\bigr{)}u_{\varepsilon}\bigr{|}\leqslant\frac{C}{\varepsilon^{4}}|u_{\varepsilon}|^{3}\,,$

it follows by Cauchy-Schwarz that

$\displaystyle\int\eta^{2}|\nabla u_{\varepsilon}-\nabla
v_{\varepsilon}|^{2}+\frac{c_{0}}{2\varepsilon^{4}}\int\eta^{2}|u_{\varepsilon}-v_{\varepsilon}|^{2}\,dx$
(2.41) $\displaystyle\quad\quad\leqslant
C\biggl{(}\frac{1}{\varepsilon^{4}}\int\eta^{2}|v_{\varepsilon}-u_{\varepsilon}|^{2}\biggr{)}^{\frac{1}{2}}\biggl{(}\frac{1}{\varepsilon^{4}}\int\eta^{2}|u_{\varepsilon}|^{6}\biggr{)}^{\nicefrac{{1}}{{2}}}+C\int|\nabla\eta|^{2}|v_{\varepsilon}-u_{\varepsilon}|^{2}\,dx$
$\displaystyle\quad\quad\leqslant\frac{c_{0}}{4\varepsilon^{4}}\int\eta^{2}|u_{\varepsilon}-v_{\varepsilon}|^{2}\,dx+\frac{C}{\varepsilon^{4}}\int\eta^{2}|u_{\varepsilon}|^{6}\,dx+C\int|\nabla\eta|^{2}|v_{\varepsilon}-u_{\varepsilon}|^{2}\,dx\,,$

so that upon rearranging, we get

$\int\eta^{2}|\nabla u_{\varepsilon}-\nabla
v_{\varepsilon}|^{2}+\frac{c_{0}}{4\varepsilon^{4}}\int\eta^{2}|u_{\varepsilon}-v_{\varepsilon}|^{2}\,dx\lesssim\frac{C}{2\varepsilon^{4}}\int\eta^{2}|u_{\varepsilon}|^{2}+C\int|\nabla\eta|^{2}|v_{\varepsilon}-u_{\varepsilon}|^{2}\,dx\,.$
(2.42)

Let $\alpha>0$ be a parameter that will be subsequently chosen. Further, we
let $\chi\in C^{\infty}_{c}(\mathbb{R}^{2}),$ be a test function that will be
chosen, which will satisfy $0\leqslant\chi\leqslant 1,$ and $\chi\equiv 0$ in
the ball $B_{4}$ of radius $4$ centered at the origin. Finally, we set $\eta$
is defined by

$\eta(x)=e^{\alpha H_{\varepsilon}}\chi\,,$

where $H_{\varepsilon}$ is given in (2.40). As $|g|\leqslant 1,$ we note that
$|u_{\varepsilon}(x)|\lesssim e^{-H_{\varepsilon}}$ (for example by examining
the representation formula using the Green’s function of the operator
$\Delta-\varepsilon^{-4}$), and that on the support of $\chi,$ we have

$|\nabla\eta|\leqslant e^{\alpha H_{\varepsilon}}\bigl{(}\alpha\chi|\nabla
H_{\varepsilon}|+|\nabla\chi|\bigr{)}\leqslant e^{\alpha
H_{\varepsilon}}\biggl{(}\alpha\frac{\chi}{\varepsilon^{2}}+|\nabla\chi|\biggr{)}=\frac{\alpha}{\varepsilon^{2}}\eta+e^{\alpha
H_{\varepsilon}}|\nabla\chi|\,.$

Inserting this in the estimate (2.42) we obtain

$\displaystyle\int\eta^{2}|\nabla u_{\varepsilon}-\nabla
v_{\varepsilon}|^{2}+\frac{c_{0}}{4\varepsilon^{4}}\eta^{2}|u_{\varepsilon}-v_{\varepsilon}|^{2}\,dx\lesssim\frac{1}{\varepsilon^{4}}\int\eta^{2}|u_{\varepsilon}|^{6}\,dx+\int\biggl{(}\frac{\alpha^{2}}{\varepsilon^{4}}\eta^{2}+e^{2\alpha
H_{\varepsilon}}|\nabla\chi|^{2}\biggr{)}|v_{\varepsilon}-u_{\varepsilon}|^{2}\,dx\,.$

Choosing $\alpha^{2}=\max(1,\frac{c_{0}}{8}),$ we obtain

$\int e^{2\alpha H_{\varepsilon}}\chi^{2}|\nabla u_{\varepsilon}-\nabla
v_{\varepsilon}|^{2}+\frac{c_{0}}{8\varepsilon^{4}}\int e^{2\alpha
H_{\varepsilon}}\chi^{2}|u_{\varepsilon}-v_{\varepsilon}|^{2}\lesssim\frac{1}{\varepsilon^{4}}\int
e^{(2\alpha-6)H_{\varepsilon}}\chi^{2}\,dx+\int e^{2\alpha
H_{\varepsilon}}|\nabla\chi|^{2}|v_{\varepsilon}-u_{\varepsilon}|^{2}\,dx\,.$

To conclude, we simply choose a sequence of dyadic annuli
$\mathcal{A}_{k}:=\\{|x|\in[2^{k-1}R_{0},2^{k+2}R_{0}]\\}$, and a
corresponding sequence of choices $\chi=\chi_{k}$, with $\chi_{k}(x)\equiv 1$
for $|x|\in[2^{k}R_{0},2^{k+1}R_{0}],$ and $\chi_{k}(x)\equiv 0$ when
$|x|\leqslant 2^{k-1}R_{0}$ or if $|x|\geqslant 2^{k+2}R_{0},$ and
$|\nabla\chi_{k}|\leqslant\frac{C}{R_{0}2^{k}}.$ Summing over
$k\in\mathbb{N},$ and buckling one last time, we find

$\int_{|x|\geqslant R_{0}}e^{2\alpha H_{\varepsilon}}\biggl{(}|\nabla
u_{\varepsilon}-\nabla
v_{\varepsilon}|^{2}+\frac{c_{0}}{16\varepsilon^{4}}|u_{\varepsilon}-v_{\varepsilon}|^{2}\biggr{)}\,dx\lesssim\frac{C}{\varepsilon^{4}}e^{-4R_{0}}+\frac{C}{R_{0}^{2}}\int_{|x|\leqslant
R_{0}}e^{2\alpha H_{\varepsilon}}|u_{\varepsilon}-v_{\varepsilon}|^{2}\,.$

Finally, since $|u_{\varepsilon}|\leqslant 1$ and $|v_{\varepsilon}|\leqslant
1,$ we conclude by the triangle inequality and multiplying through by
$\varepsilon^{4}$, that

$\int_{|x|\geqslant R_{0}}e^{2\alpha
H_{\varepsilon}}|u_{\varepsilon}-v_{\varepsilon}|^{2}\,dx\lesssim
Ce^{-4R_{0}}+C\varepsilon^{4}e^{2\alpha R_{0}}\,.$

The desired estimate follows upon optimizing in $R_{0}.$ ∎

## 3\. Two particles and nonconstant boundary conditions

In this section, we continue working with the geometry of Section 2, but
consider, instead, variable boundary conditions. In other words, we seek to
obtain an expansion, in powers of $\varepsilon$ of the energy
$F_{\varepsilon}$ for the problem (1.6), when $g_{1},g_{2}$ are nonconstant.
Without loss of generality, we assume that the functions $g_{i}:\partial
B_{i}^{\varepsilon}\to\mathbb{R}^{2}$ have uniformly and absolutely summable
Fourier developments with respect to local polar coordinates on the circles
$\partial B_{i}^{\varepsilon}.$ To be precise, parametrizing $\partial B_{1}$
by
$\\{(0,\nicefrac{{1}}{{\varepsilon^{2}}}+b)+\nicefrac{{1}}{{\varepsilon^{2}}}(\cos\theta,\sin\theta):\theta\in[0,2\pi)\\},$
we assume that

$\displaystyle
g_{1}\bigl{(}(0,\nicefrac{{1}}{{\varepsilon^{2}}}+b)+\nicefrac{{1}}{{\varepsilon^{2}}}(\cos\theta,\sin\theta)\bigr{)}=\sum_{m\in\mathbb{Z}}g_{m}^{(1)}e^{im\theta},\quad\quad\theta\in[0,2\pi),$

for some Fourier coefficients
$\\{g_{m}^{(1)}\\}_{m\in\mathbb{Z}}\in\ell^{2}(\mathbb{Z};\mathbb{C}^{2}).$
Similarly, parametrizing $\partial B_{2}$ by
$\\{(0,-\nicefrac{{1}}{{\varepsilon^{2}}}-b)+\nicefrac{{1}}{{\varepsilon^{2}}}(\cos\theta,\sin\theta):\theta\in[0,2\pi)\\},$
we assume that

$\displaystyle
g_{2}\bigl{(}(0,-\nicefrac{{1}}{{\varepsilon^{2}}}-b)+\nicefrac{{1}}{{\varepsilon^{2}}}(\cos\theta,\sin\theta)\bigr{)}=\sum_{m\in\mathbb{Z}}g_{m}^{(2)}e^{im\theta},\quad\quad\theta\in[0,2\pi),$

for some Fourier coefficients
$\\{g_{m}^{(2)}\\}_{m\in\mathbb{Z}}\in\ell^{2}(\mathbb{Z};\mathbb{C}^{2}).$

As in the case of constant boundary conditions, we perform a splitting of the
solution. To be precise, we introduce, for $j=1,2,$ the functions

$\displaystyle\Psi_{m}^{(j)}(x):=\frac{K_{m}(|x-a_{j}|)}{K_{m}(\nicefrac{{1}}{{\varepsilon^{2}}})}e^{im\theta_{j}},\quad
x\in\mathbb{R}^{2}\setminus\overline{B}_{1}\cup\overline{B}_{2}$

where $\theta_{j}\in[0,2\pi)$ denotes the polar angle with pole $a_{j}$ so
that for any $x\neq a_{j},$ we have $(x-a_{j})=|x-a_{j}|e^{i\theta_{j}}.$
Here, $K_{m}$ denotes the modified Bessel function of second kind of order
$m,$ and by Lemma B.1, the function $\Psi_{m}^{(j)}$ captures the behavior of
a single colloid. Finally, as in the case of the constant boundary conditions,
we define $Z_{m}^{(j)}\in
H^{1}\mathopen{}\mathclose{{}\left(\Omega_{\varepsilon}}\right),\ j=1,2,\
m\in\mathbb{N}$ to be the unique solutions to the problems

$\displaystyle\Delta Z_{m}^{(j)}$ $\displaystyle=Z_{m}^{(j)}\ \
\quad\quad\mbox{ in }\Omega_{\varepsilon},$ (3.1) $\displaystyle Z_{m}^{(j)}$
$\displaystyle=0\quad\ \ \ \quad\quad\mbox{ on }\partial B_{j}$ $\displaystyle
Z_{m}^{(j)}$ $\displaystyle=-\Psi^{(j)}_{m}\quad\quad\mbox{ on }\partial
B_{\sigma(j)}\,.$

where we recall that the transposition map $\sigma:\\{1,2\\}\to\\{1,2\\}$ is
defined via $\sigma(1)=2,\sigma(2)=1.$

Then, it is clear that the unique solution to (1.6) is given by the formula

$\displaystyle
U=\sum_{j=1}^{2}\sum_{m\in\mathbb{Z}}g_{m}^{(j)}(\Psi_{m}^{(j)}+Z_{m}^{(j)})$
(3.2)

We will see shortly that the infinite sum in (3.2) does indeed converge in
$H^{1}(\Omega_{\varepsilon})$ and is therefore well-defined. In order to focus
on the essential issues for the time being, let us suppose that there exists
$M_{0}\in\mathbb{N},$ such that

$g_{m}^{(j)}\equiv 0\quad\quad\mbox{ for all }|m|\geqslant M_{0},j=1,2.$ (3.3)

This makes the infinite sums in (3.2), in fact, finite. With this assumption,
in what follows we will freely interchange various integrals and sums, keeping
careful track of the dependence of errors on the tail parameter $M_{0},$ and
send $M_{0}\to\infty$ at the end of the proof of Theorem 1.3 below.

Next, let us note that identical to (2.5), in this case, too, the energy
associated to $u$ admits a splitting. To be precise, we have

######  Lemma 3.1.

Under the assumption (3.3), we have the decomposition

$\displaystyle
F_{\varepsilon}(U)=\mathfrak{Re}\mathopen{}\mathclose{{}\left[\begin{pmatrix}\ldots
g_{m}^{(1)}\ldots|\ldots
g_{m}^{(2)}\ldots\end{pmatrix}\mathcal{A}_{M_{0}}\begin{pmatrix}\vdots\\\
\overline{g_{n}^{(1)}}\vdots\\\ \---\\\
\vdots\overline{g_{n}^{(2)}}\vdots\end{pmatrix}}\right],$ (3.4)

with the matrix
$\mathcal{A}_{M_{0}}\in\mathbb{R}^{(4M_{0}+1)\times(4M_{0}+1)}$ being given by

$\begin{pmatrix}\begin{matrix}\ddots&\ddots&\ddots\\\ &&\\\
\ddots&\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(1)}+Z_{n}^{(1)})}\,d\mathcal{H}^{1}&\ddots\\\
&&\\\
\ddots&\ddots&\ddots\end{matrix}&\vline&\begin{matrix}\ddots&\ddots&\ddots\\\
&&\\\ \ddots&\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}&\ddots\\\
&&\\\ \ddots&\ddots&\ddots\end{matrix}\\\ \hline\cr\\\
\begin{matrix}\ddots&\ddots&\ddots\\\ &&\\\ \ddots&\int_{\partial
B_{2}}e^{im\theta_{2}}\overline{\frac{\partial}{\partial\nu_{2}}(\Psi_{n}^{(1)}+Z_{n}^{(1)})}\,d\mathcal{H}^{1}&\ddots\\\
&&\\\
\ddots&\ddots&\ddots\end{matrix}&\vline&\begin{matrix}\ddots&\ddots&\ddots\\\
&&\\\ \ddots&\int_{\partial
B_{2}}e^{im\theta_{2}}\overline{\frac{\partial}{\partial\nu_{2}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}&\ddots\\\
&&\\\ \ddots&\ddots&\ddots\end{matrix}\end{pmatrix}\,.$

###### Proof.

Indeed, plugging in the representation formula (3.2) for the solution $U$ into
the energy and integrating by parts, we arrive at

$\displaystyle 2F_{\varepsilon}(U)=\int_{\Omega_{\varepsilon}}|\nabla
U|^{2}+|U|^{2}\,dx=2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}U\cdot\overline{\frac{\partial
U}{\partial\nu_{1}}}\,d\mathcal{H}^{1}+\int_{\partial
B_{2}}U\cdot\overline{\frac{\partial
U}{\partial\nu_{2}}}\,d\mathcal{H}^{1}\biggr{]}.$ (3.5)

For each $j=1,2,$ on the boundary $\partial B_{j},$ we note that
$U=\sum_{m\in\mathbb{Z}}g_{m}^{(j)}e^{im\theta_{j}},$ since on $\partial
B_{\sigma(j)}$ we have $\Psi_{m}^{(j)}+Z_{m}^{(j)}=0$ by construction.
Inserting the Fourier development into (3.5), and rewriting in as a quadratic
form with the matrix being written in block form, we find

$\displaystyle 2F_{\varepsilon}(U)=2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}\biggl{(}\sum_{m\in\mathbb{Z}}g_{m}^{(1)}e^{im\theta_{1}}\biggr{)}\cdot\overline{\biggl{(}\sum_{n\in\mathbb{Z}}\sum_{j=1}^{2}g_{n}^{(j)}\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(j)}+Z_{n}^{(j)})\biggr{)}}\,d\mathcal{H}^{1}\biggr{]}$
(3.6) $\displaystyle\quad\quad+2\mathfrak{Re}\biggl{[}\int_{\partial
B_{2}}\biggl{(}\sum_{m\in\mathbb{Z}}g_{m}^{(2)}e^{im\theta_{2}}\biggr{)}\cdot\overline{\biggl{(}\sum_{n\in\mathbb{Z}}\sum_{j=1}^{2}g_{n}^{(j)}\frac{\partial}{\partial\nu_{2}}(\Psi_{n}^{(j)}+Z_{n}^{(j)})\biggr{)}}\,d\mathcal{H}^{1}\biggr{]}\,.$

Expanding and rewriting in matrix form completes the proof of the lemma. ∎

As before, our main task is to estimate the asymptotics as $\varepsilon\to
0^{+}$ of the entries of the matrix $\mathcal{A}_{M_{0}}.$ We accomplish this
in a series of Lemmas. Our first lemma is the analog of Lemma 2.2 for the
present nonconstant boundary conditions case and has a similar proof, as we
demonstrate.

######  Lemma 3.2.

The functions $Z_{n}^{(j)}$ satisfy the estimate

$\displaystyle\biggl{\|}\frac{\partial
Z_{n}^{(j)}}{\partial\nu_{j}}\biggr{\|}_{H^{-\nicefrac{{1}}{{2}}}(\partial
B_{j})}+\biggl{\|}\frac{\partial
Z_{n}^{(j)}}{\partial\nu_{\sigma(j)}}\biggr{\|}_{H^{-\nicefrac{{1}}{{2}}}(\partial
B_{\sigma(j)})}\leqslant\frac{C(1+|n|)^{\nicefrac{{1}}{{2}}}}{\sqrt{\varepsilon}}K_{n}(2b).$
(3.7)

###### Proof.

As mentioned before, the proof of this lemma proceeds similarly to that of
Lemma 2.2. Without loss of generality, fix $j=1.$
By the definition of $H^{-\nicefrac{{1}}{{2}}}(\partial\Omega_{\varepsilon})$
(see Appendix A) we have

$\displaystyle\biggl{\|}\frac{\partial
Z_{n}^{(1)}}{\partial\nu}\biggr{\|}_{H^{-\nicefrac{{1}}{{2}}}(\partial\Omega_{\varepsilon})}=\sup_{\phi\in
H^{1}(\Omega_{\varepsilon}):\|\phi\|_{H^{1}(\Omega_{\varepsilon})}\leqslant
1}\biggl{[}\int_{\Omega_{\varepsilon}}\nabla Z_{n}^{(1)}\cdot\nabla\phi+\Delta
Z_{n}^{(1)}\phi\,dx\biggr{]}$ (3.8) $\displaystyle\quad=\sup_{\phi\in
H^{1}(\Omega_{\varepsilon}):\|\phi\|_{H^{1}(\Omega_{\varepsilon})}\leqslant
1}\biggl{[}\int_{\Omega_{\varepsilon}}\nabla
Z_{n}^{(1)}\cdot\nabla\phi+Z_{n}^{(1)}\phi\,dx\biggr{]}$
$\displaystyle\quad\leqslant\|Z_{n}^{(1)}\|_{H^{1}(\Omega_{\varepsilon})}=\sqrt{2F_{\varepsilon}(Z_{n}^{(1)})},$

and therefore, the desired estimate follows by constructing a competitor to
the variational problem of minimizing the energy $F_{\varepsilon}$ subject to
the boundary conditions of $Z_{n}^{(1)}.$ Our competitor construction and
estimation of its energy proceeds as before.

_Construction of competitor $\zeta$ for $Z_{n}^{(1)}$ and estimating the
energy of the $\zeta$:_ Our competitor $\zeta\in H^{1}(\Omega_{\varepsilon})$
must be constructed satisfying the boundary conditions for $\zeta=Z_{n}^{(1)}$
on $\partial\Omega_{\varepsilon},$ so that $\zeta=0$ on $\partial B_{1}$ and
$\zeta=-\Psi_{n}^{(1)}$ on $\partial B_{2}.$ We let $\eta:(0,\infty)\to[0,1]$
be a $C^{1}$ function that satisfies $\eta(t)\equiv 1$ for
$t\in\bigl{[}\nicefrac{{1}}{{\varepsilon^{2}}},\nicefrac{{1}}{{\varepsilon^{2}}}+\nicefrac{{b}}{{2}}\bigr{]},$
$\eta(t)\equiv 0$ when $t\geqslant\nicefrac{{1}}{{\varepsilon^{2}}}+b,$ and
$|\eta^{\prime}|\leqslant\frac{2}{b},$ and set

$\displaystyle\zeta(x):=-\Psi_{n}^{(1)}(x)\eta(|x-a_{2}|).$

Then

$\displaystyle\nabla\zeta=-\eta(|x-a_{2}|)\nabla\Psi_{n}^{(1)}(x)-\Psi_{n}^{(1)}(x)\eta^{\prime}(|x-a_{2}|)\frac{x-a_{2}}{|x-a_{2}|},$

so that, pointwise, we have the bound

$\displaystyle|\nabla\zeta(x)|\leqslant|\nabla\Psi_{n}^{(1)}(x)|+\frac{2}{b}|\Psi_{n}^{(1)}(x)|,$

with support in the set
$\nicefrac{{1}}{{\varepsilon^{2}}}\leqslant|x-a_{2}|\leqslant\nicefrac{{1}}{{\varepsilon^{2}}}+b.$
Then, the energy of $\zeta$ is easily calculated:

$\displaystyle F_{\varepsilon}(\zeta)\leqslant
C\int_{\nicefrac{{1}}{{\varepsilon^{2}}}\leqslant|x-a_{2}|\leqslant\nicefrac{{1}}{{\varepsilon^{2}}}+b}\bigl{(}|\Psi_{n}^{(1)}|^{2}+|\nabla\Psi_{n}^{(1)}|^{2}\bigr{)}\,dx$
(3.9)
$\displaystyle\quad=C\biggl{|}2\mathfrak{Re}\biggl{[}-\int_{|x-a_{2}|=\nicefrac{{1}}{{\varepsilon^{2}}}}\Psi_{n}^{(1)}(x)\overline{\frac{\partial\Psi_{n}^{(1)}(x)}{\partial\nu}}\,d\mathcal{H}^{1}+\int_{|x-a_{2}|=\nicefrac{{1}}{{\varepsilon^{2}}}+b}\Psi_{n}^{(1)}(x)\overline{\frac{\partial\Psi_{n}^{(1)}(x)}{\partial\nu}}\,d\mathcal{H}^{1}\biggr{]}\biggr{|},$

where we plugged in the PDE satisfied by $\Psi_{n}^{(1)}$ and integrated by
parts as before; the signs in front of the boundary integrals reflect our
choice that the corresponding unit normals point towards $a_{2}.$ Recalling
that

$\displaystyle\Psi_{n}^{(1)}(x)=\frac{K_{n}(|x-a_{1}|)}{K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}})}e^{in\theta_{1}},$

it is clear that

$\displaystyle|\nabla\Psi_{n}^{(1)}(x)|\leqslant\frac{|K_{n}^{\prime}(|x-a_{1}|)|}{K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}})}+\frac{|n|}{|x-a_{1}|}\frac{K_{n}(|x-a_{1}|)}{K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}})},\quad\quad
x\in\mathbb{R}^{2}\setminus\overline{B}_{2}.$

_Estimating the first boundary integral in ( 3.9)._ For this term, as before,
we estimate
$|\frac{\partial\Psi_{n}^{(1)}}{\partial\nu}|\leqslant|\nabla\Psi_{n}^{(1)}|,$
and use the geometric observation that when
$|x-a_{2}|=\nicefrac{{1}}{{\varepsilon^{2}}},$ parametrizing
$x=a_{2}+\nicefrac{{1}}{{\varepsilon^{2}}}(\cos\theta,\sin\theta),\theta\in(-\pi,\pi]$
with $\theta=0$ along the vertical, we have the lower bound

$\displaystyle|x-a_{1}|\geqslant\nicefrac{{1}}{{\varepsilon^{2}}}+2b+\nicefrac{{1}}{{\varepsilon^{2}}}(1-\cos\theta).$

As $K_{n}$ and $-K_{n}^{\prime}$ are both monotone decreasing functions, and
moreover, since on the circle $|x-a_{2}|=\nicefrac{{1}}{{\varepsilon^{2}}},$
we have $\,d\mathcal{H}^{1}=\nicefrac{{1}}{{\varepsilon^{2}}}\,d\theta$ in the
above parametrization, we find

$\displaystyle\biggl{|}\int_{|x-a_{2}|=\nicefrac{{1}}{{\varepsilon^{2}}}}\Psi_{n}^{(1)}\overline{\frac{\partial\Psi_{n}^{(1)}(x)}{\partial\nu}}\,d\mathcal{H}^{1}\,d\mathcal{H}^{1}\biggr{|}$
$\displaystyle\quad\leqslant\frac{2}{\varepsilon
K_{n}^{2}(\nicefrac{{1}}{{\varepsilon^{2}}})}\Biggl{|}\int_{0}^{\pi}K_{n}\bigl{(}\nicefrac{{1}}{{\varepsilon^{2}}}+2b+\nicefrac{{1}}{{\varepsilon^{2}}}(1-\cos\theta)\bigr{)}\biggl{[}-K_{n}^{\prime}(\nicefrac{{1}}{{\varepsilon^{2}}}+2b+\nicefrac{{1}}{{\varepsilon^{2}}}(1-\cos\theta))$
$\displaystyle\quad\quad\quad\quad+\frac{|n|}{\nicefrac{{1}}{{\varepsilon^{2}}}+2b+\nicefrac{{1}}{{\varepsilon^{2}}}(1-\cos\theta)}K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}}+2b+\nicefrac{{1}}{{\varepsilon^{2}}}(1-\cos\theta){\color[rgb]{1,0,0})}\biggr{]}\,d\theta\Biggr{|}$
$\displaystyle\quad\leqslant\frac{C}{\varepsilon}(1+|n|)K_{n}^{2}(2b),$

by an easy computation similar to that in the proof of Lemma 2.2. By a similar
argument, the second term in (3.9) satisfies the bound

$\biggl{|}\int_{|x-a_{2}|=\nicefrac{{1}}{{\varepsilon^{2}}}+b}\Psi_{n}^{(1)}(x)\overline{\frac{\partial\Psi_{n}^{(1)}(x)}{\partial\nu}}\,d\mathcal{H}^{1}\biggr{|}\leqslant\frac{C}{\varepsilon}(1+|n|)K_{n}^{2}(2b).$

Putting these together with the bound in (3.9) and (3.8), the proof of the
Lemma is completed. ∎

In the next lemma we obtain the $\varepsilon\to 0^{+}$ asymptotic expansion
for the diagonal blocks in the matrix in (3.4).

######  Lemma 3.3.

For every $m,n\in\mathbb{Z},$ the $(m,n)$ entry of the diagonal block of the
quadratic form (3.4) satisfy the expansion

$\displaystyle 2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(1)}+Z_{n}^{(1)})}\,d\mathcal{H}^{1}\biggr{]}=-\frac{4\pi}{\varepsilon^{2}}\frac{K_{n}^{\prime}(\nicefrac{{1}}{{\varepsilon^{2}}})}{K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}})}\delta_{mn}+R_{mn},$
(3.10)

where

$|R_{mn}|\leqslant\frac{C}{\varepsilon}K_{m}(2b)K_{n}(2b)\sqrt{(1+|m|)(1+|n|)}.$
(3.11)

The same expansion holds for the bottom right diagonal block.

###### Proof.

Let us note that since the normal $\nu_{1}$ points towards the center $a_{1},$
it follows that in the polar coordinates about $a_{1},$ the normal derivative
$\nicefrac{{\partial}}{{\partial\nu_{1}}}=-\nicefrac{{\partial}}{{\partial
r}},$ so that, on $\partial B_{1},$

$\displaystyle\frac{\partial}{\partial\nu_{1}}\Psi_{n}^{(1)}=-\frac{K_{n}^{\prime}(\nicefrac{{1}}{{\varepsilon^{2}}})}{K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}})}e^{in\theta_{1}},$

and, we have

$\displaystyle\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}\Psi_{n}^{(1)}}\,d\mathcal{H}^{1}=-\frac{1}{\varepsilon^{2}}\frac{K_{n}^{\prime}(\nicefrac{{1}}{{\varepsilon^{2}}})}{K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}})}\int_{0}^{2\pi}e^{i(m-n)\theta_{1}}\,d\theta=-\frac{2\pi}{\varepsilon^{2}}\frac{K_{n}^{\prime}(\nicefrac{{1}}{{\varepsilon^{2}}})}{K_{n}(\nicefrac{{1}}{{\varepsilon^{2}}})}\delta_{mn}.$

As usual, the Kronecker’s delta $\delta_{mn}:=1$ if $m=n$ and $\delta_{mn}:=0$
if $m\neq n.$

We focus on the second term, i.e., on estimating

$2\mathfrak{Re}\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}Z_{n}^{(1)}}\,d\mathcal{H}^{1}.$
(3.12)

The natural idea to estimate this is to directly use Lemma 3.2; however, this
direct estimate misses the observation that $Z_{n}^{(1)}$ vanishes on
$\partial B_{1}.$ To obtain a better estimate, we use Green’s second identity
which, specialized to our setting, asserts that for any pair of suitably
smooth functions $U_{1},U_{2}$ that decay at infinity sufficiently fast
satisfy the identity

$\int_{\partial B_{1}}U_{1}\frac{\partial
U_{2}}{\partial\nu_{1}}-U_{2}\frac{\partial
U_{1}}{\partial\nu_{1}}\,d\mathcal{H}^{1}+\int_{\partial
B_{2}}U_{1}\frac{\partial U_{2}}{\partial\nu_{2}}-U_{2}\frac{\partial
U_{1}}{\partial\nu_{2}}\,d\mathcal{H}^{1}=\int_{\Omega_{\varepsilon}}U_{1}\Delta
U_{2}-U_{2}\Delta U_{1}\,dx.$

Applying this identity to $U_{1}=\Psi_{m}^{(1)}$ and
$U_{2}=\overline{Z_{n}^{(1)}},$ and subsequently to the choice
$U_{1}=\overline{\Psi_{m}^{(1)}}$, and $U_{2}={Z_{n}^{(1)}}$, adding the
results in order to obtain the real parts in (3.12), we notice that the bulk
terms on the right-hand side cancel, these choices of $U_{1}$ and $U_{2}$ are
all equal to their respective Laplacians. We are consequently only left with
boundary integrals, and we get

$\displaystyle 2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}Z_{n}^{(1)}}\,d\mathcal{H}^{1}\biggr{]}=2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}\Psi_{m}^{(1)}\overline{\frac{\partial}{\partial\nu_{1}}Z_{n}^{(1)}}\,d\mathcal{H}^{1}\biggr{]}$
(3.13) $\displaystyle\quad\quad=2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{1}}\overline{Z_{n}^{(1)}}\,d\mathcal{H}^{1}+\int_{\partial
B_{2}}\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{1}}\overline{Z_{n}^{(1)}}-\Psi_{m}^{(1)}\overline{\frac{\partial
Z_{n}^{(1)}}{\partial\nu_{1}}}\,d\mathcal{H}^{1}\biggr{]}$
$\displaystyle\quad\quad=-2\mathfrak{Re}\biggl{[}\int_{\partial
B_{2}}\overline{\Psi_{n}^{(1)}}\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{1}}+\Psi_{m}^{(1)}\overline{\frac{\partial
Z_{n}^{(1)}}{\partial\nu_{1}}}\,d\mathcal{H}^{1}\biggr{]}$

using the boundary conditions satisfied by $Z_{n}^{(1)}.$ At this point,
estimating as before and using Lemma 3.2, it is easily seen that

$\displaystyle\biggl{|}2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}Z_{n}^{(1)}}\,d\mathcal{H}^{1}\biggr{]}\biggr{|}\leqslant\frac{C(1+|n|)^{\nicefrac{{1}}{{2}}}(1+|m|)^{\nicefrac{{1}}{{2}}}K_{m}(2b)K_{n}(2b)}{\varepsilon},$

since $\|\Psi_{m}^{(1)}\|_{H^{\nicefrac{{1}}{{2}}}(\partial
B_{2})}+\|\tfrac{\partial}{\partial\nu_{1}}\Psi_{m}^{(1)}\|_{H^{-\nicefrac{{1}}{{2}}}(\partial
B_{2})}\leqslant\frac{C}{\sqrt{\varepsilon}}(1+|m|)^{\nicefrac{{1}}{{2}}}K_{m}(2b),$
and the proof of the lemma is completed.

∎

Finally, we turn to evaluating the off-diagonal blocks in (3.4). The
evaluation of these is not as straightforward since the terms involved do not
have a straightforward separation of scales. To overcome this difficulty, we
manipulate the boundary integrals that occur in the off-diagonal blocks using
integrations by parts and the PDE solved by the functions involved, and this
provides for a representation where the terms involved do have a separation of
scales. At that point we can proceed very similarly to the proof of Lemma 2.4
in the case of constant boundary conditions.

######  Lemma 3.4.

For each $m,n\in\mathbb{Z},$ the $(m,n)$ term in each of the off-diagonal
blocks satisfies the expansion

$\displaystyle\Biggl{|}2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}\biggr{]}-c_{m+n}\frac{e^{-2b}\sqrt{\pi}}{\varepsilon}\Biggr{|}\leqslant\frac{(1+|m|)(1+|n|)e^{-4b}}{\varepsilon}\,,$
(3.14)

where

$c_{m+n}:=2\Bigl{(}T_{m}(0)T_{n}(0)-U_{m-1}(0)U_{n-1}(0)\Bigr{)}=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{cc}2&(m+n)\equiv
0(\mathrm{mod}\,4)\\\ 0&(m+n)\equiv 1(\mathrm{mod}\,4)\\\ -2&(m+n)\equiv
2(\mathrm{mod}\,4)\\\ 0&(m+n)\equiv
3(\mathrm{mod}\,4)\,.\end{array}}\right.\,.$

###### Proof.

Arguing exactly like in the proof of the Lemma 3.3 using Green’s second
identity, we find that

$\displaystyle 2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}\biggr{]}$
$\displaystyle\quad=2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}\Psi_{m}^{(1)}(x)\overline{\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}\biggr{]}$
$\displaystyle\quad=2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}\frac{\partial}{\partial\nu_{1}}\Psi_{m}^{(1)}\overline{(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}$
$\displaystyle\quad\quad+\int_{\partial
B_{2}}\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{2}}\overline{(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}-\int_{\partial
B_{2}}\Psi_{m}^{(1)}\overline{\frac{\partial}{\partial\nu_{2}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}\biggr{]}.$

By the boundary conditions of $Z_{n}^{(2)},$ the first integral on $\partial
B_{1}$ vanishes, and in the second integral, the function $Z_{n}^{(2)}$
vanishes on $\partial B_{2}.$ Therefore, it follows from the preceding display
and the definition of $\Psi_{n}^{(2)}$ on $\partial B_{2}$ that

$\displaystyle 2\mathfrak{Re}\biggl{[}\int_{\partial
B_{1}}e^{im\theta_{1}}\overline{\frac{\partial}{\partial\nu_{1}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}\biggr{]}$
(3.15) $\displaystyle=2\mathfrak{Re}\biggl{[}\int_{\partial
B_{2}}\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{2}}\overline{\Psi_{n}^{(2)}}\,d\mathcal{H}^{1}\biggr{]}$
$\displaystyle\quad\quad+2\mathfrak{Re}\biggl{[}\int_{\partial
B_{2}}\Psi_{m}^{(1)}\overline{\frac{\partial}{\partial\nu_{2}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}\biggr{]}.$

Using (3.2) and arguing as before using the
$H^{\nicefrac{{1}}{{2}}}-H^{-\nicefrac{{1}}{{2}}}$ estimate, it is clear that
the last term satisfies the bound

$\displaystyle\biggl{|}2\mathfrak{Re}\biggl{[}\int_{\partial
B_{2}}\Psi_{m}^{(1)}\overline{\frac{\partial}{\partial\nu_{2}}(\Psi_{n}^{(2)}+Z_{n}^{(2)})}\,d\mathcal{H}^{1}\biggr{]}\biggr{|}\leqslant\frac{C}{\varepsilon}\sqrt{(1+|m|)(1+|n|)}K_{m}(2b)K_{n}(2b).$
(3.16)

Therefore in order to complete the proof of the lemma it remains to evaluate
the first term on the right-hand side of (3.15). Toward this end, we
parametrize $\partial B_{2}$ via
$x=\\{(0,-\nicefrac{{1}}{{\varepsilon^{2}}}-b)+\nicefrac{{1}}{{\varepsilon^{2}}}(\cos\theta_{2},\sin\theta_{2}):\theta_{2}\in[0,2\pi)\\},$
and notice that in this parametrization $\Psi_{n}^{(2)}(x)=e^{in\theta_{2}},$
and $\nu_{2}=-(\cos\theta_{2},\sin\theta_{2}).$

The main contribution then is that it remains to evaluate the first term in
(3.15). We proceed identically as in the proof of Lemma 2.4. As in that
argument, it suffices once again, to evaluate the portion of the integral on
$\partial B_{2}^{+},$ and to do this we parametrize $\partial B_{2}^{+}$ as in
the proof of that lemma (as a function of
$x\in[0,\frac{1}{\varepsilon^{2}}])$, and split the associated integral in
$[0,\frac{M}{\varepsilon}],$ and
$[\frac{M}{\varepsilon},\frac{1}{\varepsilon^{2}}].$ We note that on $\partial
B_{2}^{+},$ we have that

$\Psi_{n}^{(2)}(\theta_{2})=e^{in\theta_{2}}=(\cos n\theta_{2},\sin
n\theta_{2})=(T_{n}(\cos\theta_{2}),U_{n-1}(\cos\theta_{2})\sin\theta_{2})\,,$

with $T_{n}$ and $U_{n}$ denoting the Tchebyshev polynomials of order
$n\in\mathbb{Z}$, of the first and second kinds, respectively (see also
Appendix D); introducing these special functions permits us to express
multiple angle trigonometric functions of $\theta_{2}$ in terms of $x$ and
$y$. Indeed, since we have

$(x,y)=(0,-\frac{1}{\varepsilon^{2}}-b)+\frac{1}{\varepsilon^{2}}(\cos\theta_{2},\sin\theta_{2})\,,$

it follows by rewriting $\Psi_{n}^{(2)}$ in terms of $x,$ that

$\Psi_{n}^{(2)}(x)=\Bigl{(}T_{n}(\varepsilon^{2}x),\varepsilon^{2}\sqrt{1-x^{2}}U_{n-1}(\varepsilon^{2}x)\Bigr{)}\,.$

Before computing $\nu_{2}\cdot\nabla\Psi_{m}^{(1)},$ we record that for
$(x,y)\in\partial B_{2}^{+},$ arguing as in Step 3 of the proof of Lemma 2.4,
we find

$\displaystyle(\cos\theta_{1},\sin\theta_{1})$
$\displaystyle=\frac{(x,y)-(0,\frac{1}{\varepsilon^{2}}+b)}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}$
$\displaystyle=\frac{\Bigl{(}x,-2b-\frac{2}{\varepsilon^{2}}+\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}\Bigr{)}}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}$
$\displaystyle=\frac{\Bigl{(}\varepsilon^{2}x,-2b-2+\sqrt{1-\varepsilon^{4}x^{2}}\Bigr{)}}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}\,.$

It follows that for $(x,y)\in\partial B_{2}^{+},$ we have

$\displaystyle\Psi_{m}^{(1)}(x,y)=\frac{K_{m}\Bigl{(}\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}\Bigr{)}}{K_{m}(\frac{1}{\varepsilon^{2}})}(\cos(m\theta_{1}),\sin(m\theta_{1}))$
$\displaystyle\quad=\frac{K_{m}\Bigl{(}\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}\Bigr{)}}{K_{m}(\frac{1}{\varepsilon^{2}})}$
$\displaystyle\quad\times\Biggl{(}T_{m}\Biggl{(}\frac{\varepsilon^{2}x}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}\Biggr{)},U_{m-1}\Biggl{(}\frac{\varepsilon^{2}x}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}\Biggr{)}\frac{y-\frac{1}{\varepsilon^{2}}-b}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}\Biggr{)}\,,$

and we recall from the computations in Step 3 of the proof of Lemma 2.4 that

$\nu_{2}(x,y)=-\varepsilon^{2}(x,y+b+\frac{1}{\varepsilon^{2}})\,.$

Recall that we have identified a complex-valued function with an
$\mathbb{R}^{2}-$valued one. For the real parts, we compute

$\displaystyle\nu_{2}\cdot\nabla(\Psi_{m}^{(1)})_{1}(x,y)$
$\displaystyle=\frac{K_{m}\Bigl{(}\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}\Bigr{)}}{K_{m}(\frac{1}{\varepsilon^{2}})}T_{m}^{\prime}\Biggl{(}\frac{\varepsilon^{2}x}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}\Biggr{)}\nu_{2}\cdot\nabla\frac{\varepsilon^{2}x}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}\Bigg{|}_{y=-b-\frac{1}{\varepsilon^{2}}+\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}}$
$\displaystyle+\frac{K_{m}^{\prime}\Bigl{(}\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}\Bigr{)}}{K_{m}(\frac{1}{\varepsilon^{2}})}T_{m}\Biggl{(}\frac{\varepsilon^{2}x}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}\Biggr{)}\nu_{2}\cdot\nabla\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}\Bigg{|}_{y=-b-\frac{1}{\varepsilon^{2}}+\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}}\,.$

A tedious computation yields that for $(x,y)\in\partial B_{2}^{+},$ we have

$\displaystyle\nu_{2}\cdot\nabla\frac{\varepsilon^{2}x}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}\Bigg{|}_{y=-b-\frac{1}{\varepsilon^{2}}+\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}}$
$\displaystyle=-(\varepsilon^{2}x,\sqrt{1-\varepsilon^{4}x^{2}})\cdot\begin{pmatrix}\frac{\varepsilon^{4}\Bigl{(}\sqrt{1-\varepsilon^{4}x^{2}}-2-2b\varepsilon^{2}\Bigr{)}^{2}-\varepsilon^{6}x^{2}}{(1+\varepsilon^{4}C_{M}(x))^{\nicefrac{{3}}{{2}}}}\\\
\frac{-\varepsilon^{6}x(-2-2b\varepsilon^{2}+\sqrt{1-\varepsilon^{4}x^{2}})}{(1+\varepsilon^{4}C_{M}(x))^{\nicefrac{{3}}{{2}}}}\end{pmatrix}$
$\displaystyle=O(\varepsilon^{4})\,.$

Another computation yields that for $(x,y)\in\partial B_{2}^{+},$

$\displaystyle\nu_{2}\cdot\nabla\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}$
$\displaystyle=-(\varepsilon^{2}x,\sqrt{1-\varepsilon^{4}x^{2}})\cdot\frac{(x,\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}-\frac{2}{\varepsilon^{2}}-2b)}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}$
$\displaystyle=-(\varepsilon^{2}x,\sqrt{1-\varepsilon^{4}x^{2}})\cdot\frac{(\varepsilon^{2}x,\sqrt{1-\varepsilon^{4}x^{2}}-2-2b\varepsilon^{2})}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}$
$\displaystyle=-\frac{1-(2+2b\varepsilon^{2})\sqrt{1-\varepsilon^{4}x^{2}}}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}\approx
1\,,$

when $|x|\leqslant\frac{C}{\varepsilon}.$

We turn to computing $\frac{\partial(\Psi_{m}^{(1)})_{2}}{\partial\nu_{2}}.$
For this we observe that

$\displaystyle\nu_{2}\cdot(\nabla\Psi_{m}^{(1)})_{2}(x,y)$
$\displaystyle=\frac{1}{K_{m}(\frac{1}{\varepsilon^{2}})}\nu_{2}\cdot\nabla\Biggl{[}K_{m}\Bigl{(}\sqrt{x^{2}+\bigl{(}y-\frac{1}{\varepsilon^{2}}-b\bigr{)}^{2}}\Bigr{)}U_{m-1}\Biggl{(}\frac{\varepsilon^{2}x}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}\Biggr{)}$
$\displaystyle\quad\quad\times\frac{y-\frac{1}{\varepsilon^{2}}-b}{\sqrt{x^{2}+\Bigl{(}y-\frac{1}{\varepsilon^{2}}-b\Bigr{)}^{2}}}\Biggr{]}$
$\displaystyle\quad\quad\approx-\frac{1}{K_{m}(\frac{1}{\varepsilon^{2}})}U_{m-1}\Biggl{(}\frac{\varepsilon^{2}x}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}\Biggr{)}\frac{\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}-\frac{2}{\varepsilon^{2}}-2b}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}$
$\displaystyle\quad\quad\times
K_{m}^{\prime}\Bigl{(}\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}\Bigr{)}\frac{(x,\sqrt{\frac{1}{\varepsilon^{4}}-x^{2}}-\frac{2}{\varepsilon^{2}}-2b)}{\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}}\cdot(\varepsilon^{2}x,\sqrt{1-\varepsilon^{4}x^{2}})$
$\displaystyle\quad\quad=-\frac{1}{K_{m}(\frac{1}{\varepsilon^{2}})}U_{m-1}\Biggl{(}\frac{\varepsilon^{4}x}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}\Biggr{)}\frac{\sqrt{1-\varepsilon^{4}x^{2}}-2-2b\varepsilon^{2}}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}$
$\displaystyle\quad\quad\quad\times
K_{m}^{\prime}\Bigl{(}\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}\Bigr{)}\frac{1-(2+2b\varepsilon^{2})\sqrt{1-\varepsilon^{4}x^{2}}}{\sqrt{1+\varepsilon^{4}C_{M}(x)}}$
$\displaystyle\quad\quad\approx-
U_{m-1}(0)\frac{K_{m}^{\prime}\Bigl{(}\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}\Bigr{)}}{K_{m}(\frac{1}{\varepsilon^{2}})}\,,$

and

$(\Psi_{n}^{(2)})_{2}(x,y)=\varepsilon^{2}\sqrt{1-x^{2}}U_{n-1}(\varepsilon^{2}x)\approx\varepsilon^{2}U_{n-1}(0)\,.$

In the above, $\approx$ means that the left and right-hand sides of the
equality differ by $O(\varepsilon^{2})$ in magnitude, as can be checked from a
straightforward calculation (by noting that $\nu_{2}\approx(0,-1),$ and so
only the derivatives along the ”radial” direction contribute)

Therefore, arguing exactly as before, and recalling that

$U_{m}(0)=T_{m}(0)=\mathopen{}\mathclose{{}\left\\{\begin{aligned} 1&\quad
m\equiv 0(\mathrm{mod}\,4)\\\ 0&\quad m\equiv 1\mbox{ or
}3(\mathrm{mod}\,4)\\\ -1&\quad m\equiv
2(\mathrm{mod}\,4)\end{aligned}}\right.\,,$

we obtain

$\displaystyle 2\mathfrak{Re}\int_{\partial
B_{2}^{+}}\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{2}}\overline{\Psi_{n}^{(2)}}\,d\mathcal{H}^{1}$
$\displaystyle=2\int_{0}^{\frac{M}{\varepsilon}}\Bigl{[}\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{2}}\overline{\Psi_{n}^{(2)}}+\overline{\frac{\partial\Psi_{m}^{(1)}}{\partial\nu_{2}}}\Psi_{n}^{(2)}\Bigr{]}\frac{1}{\sqrt{1-\varepsilon^{4}x^{2}}}\,dx+O\Bigl{(}\exp\bigl{(}-\frac{c}{\varepsilon^{2}}\bigr{)}\Bigr{)}$
$\displaystyle=2\int_{0}^{\frac{M}{\varepsilon}}\Bigl{[}\frac{\partial(\Psi_{m}^{(1)})_{1}}{\partial\nu_{2}}(\Psi_{n}^{(2)})_{1}+\frac{\partial(\Psi_{m}^{(1)})_{2}}{\partial\nu_{2}}(\Psi_{n}^{(2)})_{2}\Bigr{]}\frac{1}{\sqrt{1-\varepsilon^{4}x^{2}}}\,dx+O\Bigl{(}\exp\bigl{(}-\frac{c}{\varepsilon^{2}}\bigr{)}\Bigr{)}$
$\displaystyle=2\Bigl{(}T_{m}(0)T_{n}(0)-U_{m-1}(0)U_{n-1}(0)\Bigr{)}\int_{0}^{\frac{M}{\varepsilon}}\frac{K_{m}^{\prime}\Bigl{(}\sqrt{\frac{1}{\varepsilon^{4}}+C_{M}(x)}\Bigr{)}}{K_{m}(\frac{1}{\varepsilon^{2}})}\,dx$
$\displaystyle\approx
2\Bigl{(}T_{m}(0)T_{n}(0)-U_{m-1}(0)U_{n-1}(0)\Bigr{)}\frac{e^{-2b}\sqrt{\pi}}{\varepsilon}\,,$

where we argue exactly like in the proof of Lemma 2.4. The proof of the lemma
is completed by invoking the results in Appendix D on
$c_{m+n}:=\Bigl{(}T_{m}(0)T_{n}(0)-U_{m-1}(0)U_{n-1}(0)\Bigr{)}.$

∎

With the foregoing lemmas at hand, just like in the constant boundary case,
the proof of the main theorem of this section is then immediate.

###### Proof of Theorem 1.3.

The proof is a combination of the preceding lemmas, and sending
$M_{0}\to\infty$ (see Appendix D). We note that for fixed $M_{0}<\infty,$ from
(3.4),

$\displaystyle G^{T}\mathcal{A}_{M_{0}}G$
$\displaystyle=2\sum_{m,n}c_{m+n}g_{m}^{(1)}\overline{g_{n}^{(2)}}$
$\displaystyle=\mathfrak{R}\Bigl{(}2g_{0}^{(1)}\overline{g_{0}^{(2)}}+g_{0}^{(1)}\overline{\sum_{m=0,n\neq
0}c_{n}g_{n}^{(2)}}+g_{0}^{(2)}\overline{\sum_{n=0,m\neq
0}c_{m}g_{m}^{(1)}}+\sum_{m\neq 0,n\neq
0}c_{m+n}g_{m}^{(1)}\overline{g_{n}^{(2)}}\Bigr{)}\,.$

where, since all the sums are finite, we can freely rearrange terms in the
summation, and carry out various differentiation and integration operations
term-by-term. Now we use the definition of Fourier coefficients:

$g_{m}^{(i)}=\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}e^{-im\theta}g_{i}(\theta)\,d\theta,\quad
m\in\mathbb{Z},i\in\\{1,2\\}\,.$

Inserting this in the prior expression, we obtain

$\displaystyle G^{T}\mathcal{A}_{M_{0}}G$
$\displaystyle=2\mathfrak{R}\biggl{\\{}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\Bigl{(}\sum_{m,n\in\mathbb{Z}}c_{m+n}e^{-i(n\theta-m\phi)}\Bigr{)}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi\biggr{\\}}\,.$

Each of these is a finite sum, and so we are justified in switching
integration and summation. We compute,

$\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\Bigl{(}\sum_{m,n\in\mathbb{Z}}c_{m+n}e^{-i(n\theta-m\phi)}\Bigr{)}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi$
$\displaystyle\quad=\sum_{m,n\in\mathbb{Z}}c_{m+n}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}e^{-i(n\theta-m\phi)}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi$
$\displaystyle\quad=\sum_{m,n}c_{m+n}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}e^{-i(m+n)\theta}e^{im\theta}e^{im\phi}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi$
$\displaystyle\quad=\sum_{k,m\in\mathbb{Z}}c_{k}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}e^{-ik\theta}e^{im\theta}e^{im\phi}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi$
$\displaystyle\quad=2\sum_{k,m\in\mathbb{Z}}\biggl{[}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}e^{-i4k\theta}e^{im\theta}e^{im\phi}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi$
$\displaystyle\qquad\qquad-\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}e^{-i(4k+2)\theta}e^{im\theta}e^{im\phi}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi\biggr{]}$
$\displaystyle\quad=2\sum_{m\in\mathbb{Z}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\sum_{k\in\mathbb{Z}}e^{-i4k\theta}\bigl{(}1-e^{-2i\theta}\bigr{)}e^{im\theta}e^{im\phi}g_{1}(\theta)\overline{g_{2}(\phi)}\,d\theta\,d\phi\,.$

By Lemma 3.5 below, the summation and integration in $\theta$ can be carried
out; the last line then rewrites as

$\displaystyle\frac{1}{2}\sum_{m\in\mathbb{Z}}\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\sum_{j=1}^{4}\bigl{(}1-e^{ij\pi}\bigr{)}e^{\frac{imj\pi}{2}}g_{1}\bigl{(}\frac{j\pi}{2}\bigr{)}e^{im\phi}\overline{g_{2}(\phi)}\,d\phi$
$\displaystyle\quad=g_{1}\biggl{(}\frac{\pi}{2}\biggr{)}\overline{g_{2}\biggl{(}-\frac{\pi}{2}\biggr{)}}+g_{1}\biggl{(}\frac{3\pi}{2}\biggr{)}\overline{g_{2}\biggl{(}-\frac{3\pi}{2}\biggr{)}}$
$\displaystyle\quad=2g_{1}\biggl{(}\frac{3\pi}{2}\biggr{)}\overline{g_{2}\biggl{(}-\frac{3\pi}{2}\biggr{)}}\,,$

where the last line follows by a direct computation (by, for example, writing
each of $g_{1}(\frac{\pi}{2})$ and $g_{2}(-\frac{\pi}{2})$ using Fourier
inversion, rearranging, and simplifying). As the point $\theta=\frac{3\pi}{2}$
corresponds to the bottom tip of the upper circle $\partial B_{1}$, and
$\phi=-\frac{3\pi}{2}$ corresponds to the upper tip of the lower circle
$\partial B_{2}$, the proof of the theorem is completed when boundary
conditions $g_{1}$ and $g_{2}$ have no more than the first $M_{0}$ modes in
Fourier space. Sending $M_{0}\to\infty$ completes the proof. ∎

In the course of the above argument, we used the following elementary lemma
with the choice $p=4$:

######  Lemma 3.5.

For any continuous, $2\pi-$periodic function $f:[0,2\pi)\to\mathbb{C},$ with
absolutely convergent Fourier series, and any positive integer $k$ we have

$\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{0}^{2\pi}\sum_{p\in\mathbb{Z}}e^{-ikp\theta}f(\theta)\,d\theta=\frac{1}{k}\sum_{j=1}^{k}f\biggl{(}\frac{2\pi
j}{k}\biggr{)}\,.$

###### Proof.

We start from the classical observation that under the given assumptions on
$f$,

$f(\theta)=\sum\limits_{n\in\mathbb{Z}}\exp(in\theta)\widehat{f}_{n},$ (3.17)

where the series is absolutely convergent and uniform in $\theta$. Thus, it is
straightforward to verify that

$\begin{split}\frac{1}{k}\sum\limits_{j=1}^{k}f\mathopen{}\mathclose{{}\left(\frac{2\pi
j}{k}}\right)=&\frac{1}{k}\sum\limits_{j=1}^{k}\sum\limits_{n\in\mathbb{Z}}\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
nj}{k}}\right)\widehat{f}_{n}=\sum\limits_{n\in\mathbb{Z}}\mathopen{}\mathclose{{}\left(\frac{1}{k}\sum\limits_{j=1}^{k}\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
nj}{k}}\right)}\right)\widehat{f}_{n}.\end{split}$ (3.18)

For the interior sum, we notice that if $n$ is an integer multiple of $k$, so
that $n=pk$, then $\exp\mathopen{}\mathclose{{}\left(2\pi pj}\right)=1$, and
thus
$\frac{1}{k}\sum\limits_{j=1}^{k}\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
nj}{k}}\right)=1$. However, if $n\equiv r\not\equiv 0\pmod{k}$, we have

$\begin{split}\sum\limits_{j=1}^{k}\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
nj}{k}}\right)=&\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
r}{k}}\right)\sum\limits_{j=0}^{k-1}\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
r}{k}}\right)^{j}\\\ =&\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
r}{k}}\right)\frac{\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
r}{k}}\right)^{k}-1}{\exp\mathopen{}\mathclose{{}\left(\frac{2\pi
r}{k}}\right)-1}=0.\end{split}$ (3.19)

Thus yielding

$\frac{1}{k}\sum\limits_{j=1}^{k}f\mathopen{}\mathclose{{}\left(\frac{2\pi
j}{k}}\right)=\sum\limits_{n\in\mathbb{Z}}\delta_{n\equiv
0(\text{mod}k)}\widehat{f}_{n}=\sum\limits_{p\in\mathbb{Z}}\widehat{f}_{pk}$
(3.20)

∎

## 4\. Interaction energies of multiple particles

In this subsection we demonstrate how the analysis of the present paper can be
extended to multiple particles. We will also indicate how to modify the
arguments to permit polydisperse collections of particles. As a first step
toward these generalizations, we consider unit balls $\\{B_{i}\\}_{i=1}^{N},$
with disjoint closures

$\overline{B}_{i}\cap\overline{B}_{j}=\emptyset\,.$

Denoting the center of the ball $B_{i}$ via $a_{i}^{\varepsilon},$ for any
$i,j\in\\{1,\cdots,N\\}$ we define $b_{ij}>0$ via

$b_{ij}:=\frac{|a_{i}^{\varepsilon}-a_{j}^{\varepsilon}|-2}{\varepsilon^{2}}\,.$

We consider obtaining an energy expansion to the solution $u_{\varepsilon}\in
H^{1}\bigl{(}\mathbb{R}^{2}\setminus\bigcup_{i=1}^{N}\overline{B}_{i}\bigr{)}$
of the problem

$\displaystyle\Delta u_{\varepsilon}$
$\displaystyle=\frac{1}{\varepsilon^{4}}u_{\varepsilon}\quad\mbox{ in
}\mathbb{R}^{2}\setminus\bigcup_{i=1}^{N}\overline{B}_{i}$ (4.1)
$\displaystyle u_{\varepsilon}$ $\displaystyle=g_{i}\quad\mbox{ on }\partial
B_{i}\,.$

For simplicity, we focus on the case where the boundary conditions
$g_{\varepsilon}^{i}$ are all constant; the generalization of the discussion
here to nonconstant $g_{\varepsilon}^{i}$ can then be easily carried out.

Following (2.2) we let

$\Psi_{i}(x):=\frac{K_{0}\bigl{(}\frac{|x-a_{i}^{\varepsilon}|}{\varepsilon^{2}}\bigr{)}}{K_{0}(\varepsilon^{-2})}\,,\quad
i=1,\cdots,N;x\in\mathbb{R}^{2}\setminus\bigcup_{i=1}^{N}\overline{B}_{i}\,.$

Here $a_{i}^{\varepsilon}$ is the center of the ball $B_{i}.$ We also
introduce $R_{i}\in
H^{1}\bigl{(}\mathbb{R}^{2}\setminus\bigcup_{i=1}^{N}\overline{B}_{i}\bigr{)}$
denote the unique solution to

$\displaystyle\Delta R_{i}$
$\displaystyle=\frac{1}{\varepsilon^{4}}R_{i}\quad\mbox{ in
}\mathbb{R}^{2}\setminus\bigcup_{i=1}^{N}\overline{B}_{i}$ $\displaystyle
R_{i}$ $\displaystyle=1\qquad\mbox{ on }\partial B_{i}$ $\displaystyle R_{i}$
$\displaystyle=0\qquad\mbox{ on }\partial B_{j},j\neq i\,.$

Finally, we define

$Z_{i}=R_{i}-\Psi_{i}\,.$

Then, it is clear that the unique solution to (4.1) is given by

$u_{\varepsilon}=\sum_{i=1}^{N}g_{i}R_{i}=\sum_{i=1}^{N}g_{i}(\Psi_{i}+Z_{i})\,.$

Then, by analogy with $\eqref{e.ZjPDE},$ the function $Z_{i}$ introduced here
is a solution to the linear PDE of interest, which vanishes on the $i$th ball,
and is equal to the negative of the single particle solution $\Psi_{i}$ on all
other balls. The analogy of (2.4) and (2.5) is then apparent, and we find that

$2F_{\varepsilon}(u_{\varepsilon})=\begin{pmatrix}g_{1}&\cdots&g_{N}\end{pmatrix}\begin{pmatrix}\|\Psi_{1}+Z_{1}\|_{\varepsilon}^{2}&\cdots&\langle\Psi_{1}+Z_{1},\Psi_{n}+Z_{n}\rangle_{\varepsilon}\\\
\vdots&\ddots&\vdots\\\
\langle\Psi_{n}+Z_{n},\Psi_{1}+Z_{1}\rangle_{\varepsilon}&\cdots&\|\Psi_{n}+Z_{n}\|_{\varepsilon}^{2}\end{pmatrix}\begin{pmatrix}g_{1}\\\
\vdots\\\ g_{n}\end{pmatrix}$

From this, arguing as in the two particle case, it is clear that the energy of
the minimiser is concentrated in the necks to first order in an energy
expansion: the leading order is

$\frac{2\pi}{\varepsilon^{2}}\sum_{i=1}^{N}\int_{\partial
B_{i}}|g_{i}|^{2}\,d\mathcal{H}^{1}\,,$

the next order contribution is $O(K_{0}(2b)/\varepsilon)$ arises from nearest
neighbors from the neck in between such neighbors.

The case of polydisperse particles is also similar to handle: namely, if the
particle radii vary between $[\rho_{min},\rho_{max}]$ for some
$\rho_{min},\rho_{max}=O(1)$ in $\varepsilon,$ then one simply defines
$\Psi_{i}$ and $R_{i}$ as above, accordingly.

## 5\. Numerics and comparison to nonlinear models

In this section we use numerical simulations in COMSOL, [7], to verify the
asymptotics established in the previous sections. We begin by considering the
asymptotic expansion (2.33). Recalling that the self-energy of a single
particle is given by

$\frac{2\pi}{\varepsilon^{2}}\frac{K_{1}\mathopen{}\mathclose{{}\left(\frac{1}{\varepsilon^{2}}}\right)}{K_{0}\mathopen{}\mathclose{{}\left(\frac{1}{\varepsilon^{2}}}\right)}=\frac{2\pi}{\varepsilon^{2}}\mathopen{}\mathclose{{}\left(1+\frac{1}{2}\varepsilon^{2}+o\mathopen{}\mathclose{{}\left(\varepsilon^{2}}\right)}\right)=\frac{2\pi}{\varepsilon^{2}}+\pi+o(1),$

we will replace the $O(1)$-term in (2.33) by $\pi$ and let

$\bar{\kappa}_{\varepsilon}=\frac{2\pi}{\varepsilon^{2}}\mathopen{}\mathclose{{}\left(g_{1}^{2}+g_{2}^{2}}\right)+\frac{\sqrt{\pi/2}}{\varepsilon}\mathopen{}\mathclose{{}\left[\mathrm{Li}_{\frac{1}{2}}\,\mathopen{}\mathclose{{}\left(e^{-4b}}\right)+\Theta_{4}\mathopen{}\mathclose{{}\left(e^{-2b}}\right)+e^{-4b}}\right](g_{1}^{2}+g_{2}^{2})\\\
-\frac{4\sqrt{\pi}}{\varepsilon}\mathopen{}\mathclose{{}\left[e^{-2b}+\frac{1}{\sqrt{2}}\Theta_{3}\mathopen{}\mathclose{{}\left(e^{-2b}}\right)}\right]g_{1}g_{2}+\pi\,.$
(5.1)

Fig. 5 shows the dependence of $\kappa_{\varepsilon}$ and
$\bar{\kappa}_{\varepsilon}$ on $b$ for $\varepsilon=0.15$ when $g_{1}\equiv
1$ on $\partial B_{1}$ and $g_{2}\equiv 0$ on $\partial B_{2}.$

Figure 5: Comparison between $\kappa_{\varepsilon}$ and
$\bar{\kappa}_{\varepsilon}$ when $g_{1}\equiv 1$ and $g_{2}\equiv 0$.

Fig. 6 shows the dependence of $\kappa_{\varepsilon}$ and
$\bar{\kappa}_{\varepsilon}$ on $b$ for $\varepsilon=0.15$ when $g_{1}\equiv
1$ on $\partial B_{1}$ and $g_{2}\equiv 1$ on $\partial B_{2}.$

Figure 6: Comparison between $\kappa_{\varepsilon}$ and
$\bar{\kappa}_{\varepsilon}$ when $g_{1}\equiv 1$ and $g_{2}\equiv 1$.

Fig. 7 shows the dependence of $\kappa_{\varepsilon}$ and
$\bar{\kappa}_{\varepsilon}$ on $b$ for $\varepsilon=0.15$ when $g_{1}\equiv
1$ on $\partial B_{1}$ and $g_{2}\equiv-1$ on $\partial B_{2}.$

Figure 7: Comparison between $\kappa_{\varepsilon}$ and
$\bar{\kappa}_{\varepsilon}$ when $g_{1}\equiv 1$ and $g_{2}\equiv-1$.

From Figs. 5-7 we conclude that our asymptotics are, in fact, accurate up to
$o(1).$ We also observe that for certain combinations of $g_{1}$ and $g_{2}$
the form of the $(b,\bar{\kappa}_{\varepsilon})$-dependence can be of the
Lennard-Jones-type as shown in Fig. 8.

Figure 8: $\bar{\kappa}_{\varepsilon}(b)$ when $g_{1}\equiv 1$ and
$g_{2}\equiv 0.2$.

The comparison between the energies of minimizers of the full nonlinear and
linear problems are shown in Fig. 9 for $\varepsilon=0.15$ and two different
choices of boundary data when $k(T)=2.$

Figure 9: Comparison between the energies of minimizers for the nonlinear and
linear problems when $g_{1}\equiv 1$ and $g_{2}\equiv 1$ (left) and
$g_{1}\equiv\cos{\theta}$ and $g_{2}\equiv-\cos{\theta}$ (right). Here
$u_{\varepsilon}=\operatorname*{argmin}{E_{\varepsilon}(u)}$ and
$v_{\varepsilon}=\operatorname*{argmin}{F_{\varepsilon}(u)}.$

Note that the difference between the minimum
$F_{\varepsilon}(v_{\varepsilon})$ of the quadratic energy and
$E_{\varepsilon}(u_{\varepsilon})$ of its nonlinear counterpart is roughly a
constant, hence the interaction forces between the two particles in the
nonlinear and linear regimes are approximately the same.

Next, we consider a system of three particles with $g_{1}=g_{2}=g_{3}\equiv 1$
on $\partial B_{1},$ $\partial B_{2}$ and $\partial B_{3},$ respectively. To
this end, we denote by

$\kappa_{\varepsilon}^{0}:=\frac{3\pi}{2}\mathopen{}\mathclose{{}\left(\frac{4}{\varepsilon^{2}}+1}\right)$

the self-energy of three particles and

$\kappa_{\varepsilon}^{1}(b):=\frac{\sqrt{2\pi}}{\varepsilon}\mathopen{}\mathclose{{}\left[\mathrm{Li}_{\frac{1}{2}}\,\mathopen{}\mathclose{{}\left(e^{-4b}}\right)+\Theta_{4}\mathopen{}\mathclose{{}\left(e^{-2b}}\right)+e^{-4}}\right]\\\
-\frac{4\sqrt{\pi}}{\varepsilon}\mathopen{}\mathclose{{}\left[e^{-2b}+\frac{1}{\sqrt{2}}\Theta_{3}\mathopen{}\mathclose{{}\left(e^{-2b}}\right)}\right]$

the interaction energy of a single neck between the two particles on the
distance $2\varepsilon^{2}b$ from each other (cf. (5.1)).

In the first numerical experiment, we assume that three particles are
positioned at the vertices of a equilateral triangle, where the distance
between the centers of any pair of particles is $2+2b\varepsilon^{2}.$
Assuming that the interactions are restricted to the necks, the minimum energy
of this configuration should be

$\kappa_{\varepsilon}\sim\bar{\kappa}_{\varepsilon}:=\kappa_{\varepsilon}^{0}+3\kappa_{\varepsilon}^{1}.$

Fig. 10 demonstrates that this is indeed the case as the graphs of
$\kappa_{\varepsilon}$ and $\bar{\kappa}_{\varepsilon}$ as functions of $b$
are essentially indistinguishable.

Figure 10: Comparison between $\kappa_{\varepsilon}$ and
$\bar{\kappa}_{\varepsilon}$ (left) for a three-particle triangular
configuration depicted on the right.

In Fig. 11,

Figure 11: The interaction energy between the particles forming the
configuration shown on the bottom as a function of $\theta$ (top).

we consider a configuration of three particles where the distances between
$B_{1}$ and $B_{2}$ and $B_{1}$ and $B_{3}$ are fixed and equal to
$2\varepsilon^{2},$ while the angle $\theta$ ranges from $0$ to $60$ degrees.
We plot the ratio between the interaction energy between the three particles
and the energy of a single neck. As it can be seen from Fig. 11, when the
angle is less than $60$ degrees, then there are exactly two necks and, indeed,
the interaction energy is equal to exactly two neck energies. When the angle
is equal to $60$ degrees, then the third neck forms and the interaction energy
is equal to the three energies of a single neck. We conjecture that for all
boundary conditions on the surfaces of particles, the energy of interaction
between the particles is concentrated in the necks.

## 6\. Monte Carlo simulation of the pairwise energy

In this section we will consider a many-body system of identical particles
satisfying given canonical degree $d$ boundary data
$g_{k}(\theta)=e^{i(d\theta-(d-1)\omega_{k})}$ and which we assume are free to
move and rotate, so that each particle has degrees of freedom corresponding to
its centre of mass and the angle $\omega_{k}$. Following our discussion in the
preceding section, we assume that there are only pairwise interactions in this
system. Assuming that the distances between particles are large enough and
plugging in special choices of $g_{1}(\theta)=e^{i(d\theta-(d-1)\omega_{1})}$
and $g_{2}(\theta)=e^{i(d\theta-(d-1)\omega_{2})}$ in the statement of Theorem
1.3, we obtain that the interaction energy between a pair of particles of
radius $1$ at distance $2b\varepsilon^{2}$ is given by

$V_{12}=-(-1)^{d}\frac{2e^{-2b}\sqrt{\pi}}{\varepsilon}\cos((d-1)(\omega_{1}-\omega_{2}))\,.$
(6.1)

We note a qualitative difference in behaviour depending on the parity of $d$.
When $d$ is even, the energy is minimised at parallel configurations, with
particles at relative angle of $0$. Thus, we expect it to be favourable for
particles to be closely packed with similar orientations. If $d$ is odd,
however, then the energy is minimised at anti-parallel configurations, where
particles are at a relative angle of $\frac{\pi}{d-1}$, modulo
$\frac{2\pi}{d-1}$. As the interactions are short-range, we expect only
interactions with nearest neighbours to be significant. Heuristically, it
seems clear that configurations of square-like lattices with second-nearest-
neighbours having the same orientation, whilst nearest neighbours are at a
relative angle of $\frac{\pi}{d-1}$, should be relatively stable.

We will consider the pairwise interaction energy for particles with
orientations $\omega_{1},\omega_{2}$ and centres of mass separated by $r$ as
given by

$V(r,\omega_{1},\omega_{2})=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{l
l}(-1)^{d+1}\exp\mathopen{}\mathclose{{}\left(-\frac{|r|-2}{\epsilon^{2}}}\right)\cos((d-1)(\omega_{1}-\omega_{2}))&|\widehat{r}|>2\\\
+\infty&|\widehat{r}|\leq 2\end{array}}\right.$ (6.2)

Multiplicative factors that do not affect minimisers of the energy are
neglected for simplicity. The infinite energy for $|\widehat{r}|\leq 2$
corresponds to the the particles being unable to interpenetrate. Of course,
the pairwise interaction in (6.1) corresponds to an asymptotic limit, and thus
we are required to introduce an appropriate length-scale for the simulation,
corresponding to the choice of $\epsilon$, which we take to be
$\sqrt{\frac{2}{5}}$, as we found in preliminary studies that even marginally
smaller values of $\epsilon$ lead to interactions too weak to produce any
noticable structure. The total pairwise energy is then given by

$\frac{1}{2}\sum\limits_{i\neq j}V(x_{i}-x_{j},\omega_{i},\omega_{j}).$

We employ a simulated annealing algorithm, with the transition probabilities
taken from the corresponding Gibbs’ distribution of the system, that is, those
of a Metropolis-Hastings algorithm. At each iteration step, we perturb the
centre of mass of the particle with index $i$ according to a Gaussian
distribution with mean $0$ and standard deviation given by
$\min(\max(0.025,\delta_{i}),0.5)$, where $\delta_{i}=\min\limits_{j\neq
i}(|x_{i}-x_{j}|-2)$ is the minimal contact distance to another particle. We
perturb the angle according to a normal distribution with mean 0 and standard
deviation $\frac{2\pi}{50}$. For each temperature, we trial each particle
once, and linearly decrease the temperature from $\frac{1}{4}$ to $0$ over
25000 iterations. We take a total of 256 particles, initialised as a
perturbation of a square lattice with nearest-neighbour separation of $2.2$,
and orientations taken according to a uniform distribution. Finally, due to
the short-range nature of the interactions, at higher temperatures it is easy
for particles to drift large distances, at which point their behaviour becomes
a random walk and ceases to effectively interact with the rest of the system,
so we impose that particles cannot leave a box of size $46\times 46$.

Whilst simulated annealing is generally used to find global minimisers, as the
particles are very weakly interacting, we expect a relatively flat energy
landscape that permits large fluctuations away from the global minimiser.

In Figure 1 we present the results of the simulations. The disks represent the
individual particles, and the lines within them represent their orientation,
and are illustrated such as to be consistent with the axes of symmetry of the
particles. Furthermore, we colour the particles according to their angle
modulo $|\frac{2\pi}{d-1}|$, on an RGB colour-scale, corresponding to the
symmetry of their boundary condition.

In Subfigure 1(b) we have odd-degree boundary conditions, and thus by the
interaction energy (6.1), we expect to have an anti-parallel configuration,
where neighbours in close contact are rotated, but second-nearest neighbours
have the same orientation, and this is observed. Even though a square lattice
can be expected via a heuristic argument, we observe that the distribution of
centres of mass is relatively amorphous, with some mild amount of short-range
correlation. We observe several chain-like structures, with anti-parallel
alignment with nearest neighbours. Due to the short-range nature of the
interaction, these are expected to be locally stable, whilst denser
configurations in a square-like lattice would have lower energy. This is
similar to the situation in Subfigure 1(c), which also exhibits clearly
visible well-aligned domains with anti-parallel configurations within.
Notably, we observe that the degree 5 configuration is more amorphous with
more chain-like structures, which we aim to explain via a heuristic argument.
In the degree 5 case, we need to consider relative angles modulo
$\frac{\pi}{2}$, whilst in the degree 3 case we consider relative angles
modulo $\pi$. This smaller range of angles would suggest a higher sensitivity
to small perturbations in the orientation, making it more difficult for the
particles to align into their optimum states and leading to amorphousness.

In Subfigure 1(a), our pairwise energy favours nearest neighbours having the
same orientation, and many contacts with neighbours, which is the observed
behaviour. Although particles are generally well-aligned with their
neighbours, we observe a kind of polycrystalline structure with clearly
identifiable domains. These are expected to be locally stable, as reorienting
a single grain would require simultaneous reorientation of many particles.

To demonstrate more clearly the local parallel and anti-parallel
configurations, we include histograms in Figure 12 of the relative
orientations of particles with their nearest neighbours and second-nearest
neighbours below, taken modulo $\frac{2\pi}{|d-1|}$. Explicitly, we say that
two particles are nearest neighbours if the separation of their centres of
mass is less than $2.05$, and that two particles are second-nearest neighbours
if they are distinct and share a nearest neighbour. We observe a clear
tendency for nearest neighbours to be either parallel or anti-parallel
according to the parity of the degree, and parallel alignment of next-nearest
neighbours in all cases.

(a) Degree 2: Relative angle of nearest neighbours

(b) Degree 3: Relative angle of nearest neighbours

(c) Degree 5: Relative angle of nearest neighbours

(d) Degree 2: Relative angle of next-nearest neighbours

(e) Degree 3: Relative angle of next-nearest neighbours

(f) Degree 5: Relative angle of next-nearest neighbours

Figure 12: Distributions of relative angles of nearest- and next-nearest
neigbours

It is a similarly straightforward exercise to evaluate the angular component
of the interaction energy between particles of distinct degrees. If we
consider two particles, whose centres of mass are at relative angle $\alpha$,
of degrees $d_{1},d_{2}$ and with orientations described by angles
$\omega_{1},\omega_{2}$, then we have that

$\begin{split}\mathrm{Re\,}(g_{1}(p)\bar{g_{2}}(q))=&\mathrm{Re\,}(e^{i(d_{1}\alpha-(d_{1}-1)\omega_{1})}e^{-i(d_{2}(\alpha-\pi)-(d_{2}-1)\omega_{2})})\\\
=&\cos((d_{1}\alpha-(d_{1}-1)\omega_{1})-(d_{2}(\alpha-\pi)-(d_{2}-1)\omega_{2}))\\\
=&(-1)^{d_{2}}\cos((d_{1}-d_{2})\alpha+(d_{2}-1)\omega_{2}-(d_{1}-1)\omega_{1}),\end{split}$
(6.3)

where $p$ and $q$ correspond to the closest points on the surface of each
respective particle to the other. We remark that unlike the case where both
degrees are equal, this depends on the relative position of the particles via
$\alpha$, and not just the orientations $\omega_{1},\omega_{2}$.

We consider a mixed system of degree 1 and degree 3 particles. As seen before,
we have that degree 3 particles prefer an anti-parallel alignment. Degree 1
particles are purely repulsive, and due to their rotational symmetry, there is
no orientational dependence. For the interactions between degree 1 and degree
3 particles, taking particle 1 to be of degree 3 and particle 2 to be of
degree 1, the angular component of the interaction energy is
$-\cos(2(\alpha-\omega))$. In particular, their optimal configuration is to
have the degree 1 particle at either of the two poles of the degree 3 particle
where the director is perpendicular to the surface. We employ a simulated
annealing algorithm with the same experimental setup as the previous
experiments to obtain the results in Figure 13. As before, we colour the
degree 3 particles according to their angle, modulo $\pi$, with the
illustrated diameter spanning the two points where the boundary data is
perpendicular to the surface. The degree 1 particles are rotationally
symmetric and thus coloured in white.

(a) Configuration at the end of the simulation

(b) Relative angle of the boundary data at the contact point for pairs of
degree 1 and degree 3 particles.

(c) Relative angle of the boundary data at the contact point for pairs of
degree 3 particles.

Figure 13: Results for a mixed system of degree 1 and degree 3 particles.

In Subfigure 13(a), we see the configuration at the end of the simulation. By
eye, we observe qualitatively the expected behaviour of neighbours, where
degree 1 particles are separated due to repulsive interactions, neighbouring
degree 3 particles tend to be at near-right-angles to each other, and degree 3
and degree 1 particles are roughly aligned along the illustrated diameter,
whose end-points correspond to the regions of the surface with perpendicular
director. Nonetheless, we observe that the particles are not so well-aligned
as in the pure-state case. In particular, we see many triangles consisting of
two degree 3 particles and one degree 1 particle, and geometrically such a
triangle cannot be pairwise-minimising for the energy. We demonstrate the
local orientational ordering graphically by considering the relative angles of
the boundary director at the contact point of nearest neighbours in Subfigure
13(b) for pairs of degree 1 and degree 3 particles, and in Subfigure 13(c) for
pairs of degree 3 particles, taken modulo $\pi$ in each case. We observe a
central tendency at zero, but greater variation than in the case of pure
systems.

## 7\. Acknowledgements

The authors acknowledge the hospitality of HIM: Hausdorff Center for
Mathematics where some of the research on this project was conducted during
the Trimester Program Mathematics for Complex Materials funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany
Excellence Strategy – EXC-2047/1 – 390685813. DG was supported in part by the
NSF grant DMS-2106551. R.V. was partially supported by the Simons Foundation
(Award # 733694) and an AMS-Simons travel award. He also acknowledges the
hospitality provided by the Department of Mathematics at the University of
Akron when R.V. visited D.G. to complete parts of this project. A.Z. has been
partially supported by the Basque Government through the BERC 2022-2025
program and by the Spanish State Research Agency through BCAM Severo Ochoa
excellence accreditation SEV-2017-0718 and through project
PID2020-114189RB-I00 funded by Agencia Estatal de Investigacion (PID2020-
114189RB-I00 / AEI / 10.13039/501100011033). A.Z. was also partially supported
by a grant of the Ministry of Research, Innovation and Digitization, CNCS -
UEFISCDI, project number PN-III-P4-PCE-2021-0921, within PNCDI III.

Appendix

##  Appendix A Sobolev spaces and trace theory

There are various definitions for the norms of trace spaces of functions in
$H^{1}(\Omega)$, which are equivalent for sufficiently regular $\Omega$ [8].