Desingularization of 2D elliptic free-boundary problem with non-autonomous nonlinearity

Abstract

In this paper, we consider the existence and limiting behaviour of solutions to a semilinear elliptic equation arising from confined plasma problem in dimension two

\[ \begin{cases} -\Delta u=\lambda k(x)f(u) & \text{in}\ D,\\ u= c & \displaystyle\text{on}\ \partial D,\\ \displaystyle - \int_{\partial D} \frac{\partial u}{\partial \nu}\,{\rm d}s=I, \end{cases} \]

where $D\subseteq \mathbb {R}^2$ is a smooth bounded domain, $\nu$ is the outward unit normal to the boundary $\partial D$, $\lambda$ and $I$ are given constants and $c$ is an unknown constant. Under some assumptions on $f$ and $k$, we prove that there exists a family of solutions concentrating near strict local minimum points of $\Gamma (x)=({1}/{2})h(x,\,x)- ({1}/{8\pi })\ln k(x)$ as $\lambda \to +\infty$. Here $h(x,\,x)$ is the Robin function of $-\Delta$ in $D$. The prescribed functions $f$ and $k$ can be very general. The result is proved by regarding $k$ as a $measure$ and using the vorticity method, that is, solving a maximization problem for vorticity and analysing the asymptotic behaviour of maximizers. Existence of solutions concentrating near several points is also obtained.

1. Introduction and main results

In this paper, we consider the 2D plasma problem

(1.1)\begin{equation} \begin{cases} -\Delta u=\lambda k(x)f(u) & \text{in}\ D,\\ u= c & \displaystyle\text{on}\ \partial D,\\ \displaystyle - \int_{\partial D} \frac{\partial u}{\partial \nu}\,{\rm d}s=I, \end{cases} \end{equation}

where $D\subseteq \mathbb {R}^2$ is a simply-connected bounded domain with smooth boundary, $\nu$ is the outward unit normal to the boundary $\partial D$, $\lambda \in \mathbb {R}^+$ and $I$ are given constants and $c$ is an unknown constant. The non-autonomous term $k$ and the nonlinearity $f$ are two prescribed functions. In the following, we always assume $f >0$ on $(0,\,+\infty )$ and $f \equiv 0$ on $(-\infty,\,0)$.

The problem (1.1) arises from a model describing a simplified version of equilibrium of a plasma in a cavity (also called the ‘Tokamak machine’), see [13, 35, 36]. Let $u_+=u$ if $u\geq 0$ and $u_+=0$ if $u<0$, $u_-=u_+-u$. The equilibrium of a plasma confined in a toroidal cavity is governed by the following equations (see the classical paper [35])

(1.2)\begin{equation} \begin{cases} \mathcal{L}u=\lambda k u_- & \text{in}\ D,\\ u=\gamma & \text{on}\ \partial D,\\ \displaystyle\int_{\partial D}\frac{1}{x_1}\frac{\partial u}{\partial \nu}\,{\rm d}s=I, \end{cases} \end{equation}

where $D\subseteq \mathbf {H}=\{x=(x_1,\,x_2)\mid x_1>0\}$ is a bounded domain, $\mathcal {L}=\sum _{i=1}^2({\partial }/{\partial x_i}) (({1}/{x_1})({\partial }/{\partial x_i}))$, $\nu$ is the outward unit normal to $\partial D$ at $x$, $I$ is a given positive constant and $\gamma$ is an unknown constant. The non-autonomous term $k$ satisfies

\[ 0< k_0\leq k(x)\leq k_1<\infty,\quad x\in \overline{D}. \]

The sets $D_p=\{x\in D\mid u(x)<0\}$ and $D_v=\{x\in D\mid u(x)>0\}$ are called the plasma set and vacuum set, respectively. $\partial D_p$ is called the free boundary. Indeed, from a physical point of view, the nonlinear term $\lambda k u_-$ in (1.2) can be replaced by $\lambda k f(u)$ for general $f$, which is called the $constitutive~function$ of the plasma, see Appendix in [35].

If we simplify the elliptic operator $\mathcal {L}$ in (1.2) to $\Delta$, one gets a simplified model

(1.3)\begin{equation} \begin{cases} \Delta u=\lambda k u_-\ \ & \text{in}\ D,\\ u=\gamma & \text{on}\ \partial D,\\ \displaystyle \int_{\partial D} \frac{\partial u}{\partial \nu}\,{\rm d}s=I. \end{cases} \end{equation}

Note that for a solution $u$ of (1.3), $v=-u$ satisfies (1.1) with $f(t)=t_+$ and $c=-\gamma$.

Existence of solutions to (1.3) and the general problem (1.1) is studied in many references, see [3, 13, 22, 28, 31, 33, 35] and reference therein. In the case $k\equiv 1$, it is well-known in [36] that (1.3) has a solution if and only if $\lambda >0$ and

\[ \gamma<0,\ \text{if}\ \lambda<\lambda_1;\quad \gamma=0,\ \text{if} \lambda=\lambda_1;\quad \gamma>0,\ \text{if}\ \lambda>\lambda_1, \]

where $\lambda _1$ is the first eigenvalue of $-\Delta$ in $D$ with Dirichlet boundary condition. Moreover, if $\lambda \in (0,\,\lambda _2)$, where $\lambda _2$ is the second eigenvalue of $-\Delta$ in $D$, the solution of (1.3) is unique, see also [7, 29]. For the asymptotic behaviour of solutions to (1.3) as $\lambda$ tending to infinity, Caffarelli-Friedman [10] first proved the non-uniqueness of the solutions to (1.3) and showed that the free boundary $\partial D_p$ is approximately a circle as $\lambda \to +\infty.$ Let $G(x,\,y)$ be the Green function of $-\Delta$ in $D$ with Dirichlet boundary condition. Then $G(x,\,y)$ has the decomposition

\[ G(x,y)={-}\frac{1}{2\pi}\ln|x-y|-h(x,y), \]

where $-({1}/{2\pi })\ln |x-y|$ is the fundamental solution of $-\Delta$ and the regular part $h(x,\,y)\in C^{\infty }(D\times D)$. It is proved in [10] that if there exists $O\subseteq D$ with $\min _{x\in \partial O}h(x,\, x) > \min _{x\in O}h(x,\, x)$, then (1.3) has a solution $u^\lambda$ for every $\lambda$ sufficiently large, and the corresponding plasma set $D_p$ shrinks to a point $x^*\in O$ with $h(x^*,\,x^*) = \min _{x\in O}h(x,\, x)$ as $\lambda \to +\infty$. This result has been extended to solutions to (1.3) whose plasma set has several components. Under the assumption that the homology of $\Omega$ is nontrivial, [13] proved that for every $l\geq 1$ and $\lambda$ sufficiently large, (1.3) has a solution whose plasma set $D_p$ consists of $l$ components and concentrates near critical points of the Kirchhoff–Routh Hamiltonian as $\lambda \to +\infty$ by using the Lyapunov–Schmidt reduction method. The Kirchhoff–Routh Hamiltonian $\mathcal {H}_{l}$ is defined by (see [23, 27, 30])

(1.4)\begin{equation} \mathcal{H}_{l}(x_1,\ldots,x_{l})={-}\frac{1}{2} \sum_{1\leq i\neq j\leq l}d_id_jG(x_i,x_j)+\frac{1}{2}\sum_{i=1}^{l}d_i^2h(x_i,x_i), \end{equation}

where $(x_1,\,\ldots,\,x_l)\in D^{(l)}:=\underbrace {D\times D\times \cdots \times D}_{l}$ satisfies $x_i\neq x_j$ for $i\neq j$, and $d_1,\,\ldots,\,d_l$ are $l$ prescribed constants. When $D$ has non-trivial topology, solutions of (1.1) with $f(t)=t^p_+ (p>1)$ whose plasma region shrinks down around finitely many different points have been constructed in [28]. See [4, 5, 8] for more results. Note that the plasma problem (1.1) with $k(x)\equiv 1$ also corresponds to the vorticity formulation of 2D steady incompressible Euler equations. In [12], by using the non-degeneracy of solutions to

\[ -\Delta u=u^p_+,\quad \text{in}\ \mathbb{R}^2 \]

and the Lyapunov–Schmidt finite-dimensional reduction method, [12] proved the existence of solutions of (1.1) with $f(t)=t^p_+ (p>1)$ concentrating near isolated non-degenerate critical points of $\mathcal {H}_{l}$ for $\lambda$ sufficiently large. [14] further proved the existence and asymptotic behaviour of concentrated solutions of (1.1) with $f(t)=t^p_+ (p=0)$ for $\lambda$ sufficiently large by using Lyapunov–Schmidt reduction method. Compared to [12], results in [14] require more delicate estimates since the nonlinearity in [14] is not as smooth as it is in [12]. For more results, see, e.g., [11, 15–18, 34, 37].

When $k$ is a function rather than a constant, many references also considered the existence and asymptotic behaviour of solutions to (1.1), see [19, 24, 25, 32, 38] and reference therein. [35] first obtained the existence of solutions of (1.3) by considering minimization of a certain variational problem. For $N\geq 3,$ Shibata [32] considered the following equations

(1.5)\begin{equation} \begin{cases} -\varepsilon^2\Delta u=k(x)f(u-1), \ u>0, & x\in D,\\ u=0, & x\in \partial D, \end{cases} \end{equation}

where $D\subseteq \mathbb {R}^N$, $\varepsilon >0$ is small and $k(x)$ is a positive function in $\overline {D}$. Under the assumption that $f (t) = t_+^p$ for $p\in (1,\,({N+2}/{N-2}))$, the author proved that (1.5) has a least energy solution concentrating near global maximum points of $k$ as $\varepsilon \to 0^+$. Here the concentration means that the plasma set $\{x\in D\mid u_\varepsilon (x)>1\}$ shrinks to some points as $\varepsilon \to 0$. This result has been extended to solutions to (1.5) with general nonlinearities concentrating near several boundary points, see [24]. As for the plasma problem (1.1) for $N\geq 3$, solutions whose plasma region shrinks down around finitely many distinct points were constructed in [25]. It is worth mentioning that, both in [32] and in [24, 25] the total vorticity vanishes rather than tends to a non-zero constant as $\varepsilon \to 0$, that is,

\[ \int_{\partial D} \frac{\partial u_\varepsilon}{\partial \nu}\,{\rm d}s\to 0\quad \text{as}\ \varepsilon\to0. \]

For $N=2,\,$ by considering Liouville-type equations

(1.6)\begin{equation} \begin{cases} -\Delta u=\varepsilon^2 K(x)e^u, & \text{in}\ D,\\ u=0 & \text{on}\ \partial D, \end{cases} \end{equation}

del Pino et al. [19] proved that solutions of (1.6) have the blow up-concentration phenomenon, i.e., there exist solutions of (1.6) concentrating near small neighbourhoods of critical points $\{x_{1,0},\,\ldots,\, x_{n,0}\}$ of the function

(1.7)\begin{equation} -\sum_{i=1}^n2\ln K(x_i)+8\pi h(x_i,x_i)-\sum_{j\neq i}8\pi G(x_i,x_j) \end{equation}

as $\varepsilon \to 0$, and the total vorticity of solutions around each $x_{i,0}$ tends to a non-zero constant as $\varepsilon \to 0$. Note that (1.7) is different from the Kirchhoff–Routh Hamiltonian (1.4) since the presence of $K$. Note also that (1.6) coincides with (1.1) by letting $\lambda =\varepsilon ^2$ and $f(t)=e^t$. A natural question is, whether there exist solutions $u^\lambda$ to (1.1) with general profile function $f$, such that the corresponding ‘plasma set’ $\{x\in D\mid u^\lambda (x)>0\}$ concentrates near several points with diameter tending to 0 as $\lambda \to +\infty$?

In this paper, we will construct solutions to (1.1) concentrating near some prescribed points with a large class of sub-exponential nonlinearities $f$. The nonlinearity $f$ can either be continuous (e.g., $f(t)=t^p_+$ for some $p\in (0,\,+\infty )$) or have a jump (e.g., $f$ being a Heaviside function), see theorems 1.1 and 1.3. We prove that for any $x_0$ being strict local minimizers of $\Gamma (\cdot )$ defined by (1.8), there exist solutions $u^\lambda$ of (1.1), whose ‘plasma set’ $\{x\in D\mid u^\lambda (x)>0\}$ concentrates near $x_0$ as $\lambda \to +\infty$ and total vorticity tends to a non-zero prescribed constant $I$ as $\lambda \to +\infty$. The idea is to regard the non-autonomous term $k$ as a $measure$ and to use the Arnold's variational method developed by [1, 2, 37]. Note that in [12, 14, 28], the constitutive function is $f(t)=t^p_+$ for $p\geq 0.$ Thus compared to the classical results, we can construct concentrated solutions to plasma problem (1.1) with very general nonlinearity.

Before stating our results, let us first introduce some notations: for every Lebesgue-measurable set $A,\, B\subset D$, $\overline {A}$ denotes the closure of $A$ and $|A|$ denotes the two-dimensional Lebesque measure of $A$, except when stated otherwise; $dist(A,\,B)=\inf _{x\in A, y\in B}|x-y|$ denotes the distance between $A$ and $B$; $B_{r}(y)$ denotes the open ball of radius $r$ centred at $y$; $\chi _{A}$ denotes the characteristic function of $A\in D$, namely $\chi _{A}(x)=1$ for $x\in A$ and $\chi _{A}(x)=0$ for $x\notin A$; $O(1)$ and $o(1)$ denote some quantities which remain bounded and go to zero as $\lambda$ tends to infinity, respectively; for any function $h$, $supp(h)$ denotes the support set of $h$.

Let $k(x):D\to \mathbb {R}$ be a continuous function satisfying

• ($\mathcal {K}$). There exist two constants $k_0,\, k_1$ such that

  \[ 0< k_0\leq k(x)\leq k_1<{+}\infty, \quad \forall x\in D. \]

Let us introduce an auxiliary function which will play a crucial role in the study of the existence of solutions to (1.1). Define

(1.8)\begin{equation} \Gamma(x)=\frac{1}{2}h(x,x)-\frac{1}{8\pi}\ln k(x), \end{equation}

where $h(x,\,y)$ is the regular part of the Green's function $G$, and $k$ satisfies the assumption ($\mathcal {K}$). Clearly $\Gamma (\cdot )$ is well-defined in $D.$

Our first result concerns the existence and asymptotic behaviour of solutions to (1.1) with the nonlinearity being the Heaviside function, i.e., $f(t)=\chi _{\{t>0\}}$.

Theorem 1.1 Suppose that $k(\cdot )$ satisfies $(\mathcal {K})$ and $f(t)=\chi _{\{t>0\}}$. Let $x_0$ be a strict local minimizer of $\Gamma$. Then there exists $\lambda _0>0$, such that for any $\lambda \in (\lambda _0,\,+\infty )$, (1.1) has a weak solution pair $(u^\lambda,\, c^\lambda )$ which satisfies the following properties:

1. (1) the diameter of the plasma set $\{x\in D\mid u^\lambda (x)>0\}$ is of the order $O(\lambda ^{-({1}/{2})})$ as $\lambda \to +\infty$.

2. (2) For any $x\in \{x\in D\mid u^\lambda (x)>0\}$, $x$ tends to $x_0$ as $\lambda \to +\infty$.

3. (3) For $\lambda$ sufficiently large, $\{x\in D\mid u^\lambda (x)=0\}$ is a $C^1$ curve and converges to a circle as $\lambda \to +\infty$.

4. (4) There holds

   (1.9)\begin{align} c^\lambda={-}\frac{I}{4\pi}\ln\lambda- \frac{k(x_0)}{2\pi}\int_{B_{\sqrt{{I}/{\pi k(x_0)}}}(0)}\ln\frac{1}{|x^*-y|}\,{\rm d}y+Ih(x_0,x_0)+o(1), \end{align}
   where $x^*$ is any point of $\partial B_{B_{\sqrt {{I}/{\pi k(x_0)}}}(0)}$.

Remark 1.2 We give an example to show the existence of $x_0$. By (1.8), $\Gamma (x)=({1}/{2})h(x,\,x)-({1}/{8\pi })\ln k(x)$. Since $\lim _{x\to \partial D}h(x,\,x)=+\infty$, by assumption $(\mathcal {K})$ one can get the existence of minimum points $x_{0,1}\in D$ satisfying $\Gamma (x_{0,1})=\min _{x\in D}\Gamma (x)$. Thus from theorem 1.1, there exists a family of solutions $u^\lambda$ concentrating near minimizers of $\Gamma$. Note that the limiting location of the plasma set $\{x\in D\mid u^\lambda (x)>0\}$ in theorem 1.1 coincides with that in [19] since when choosing $n=1$ in (1.7), (1.7) is equal to $16\pi$ times $\Gamma$.

When $f(t)$ is a continuous function satisfying some growth conditions, one can also get solutions to (1.1) concentrating near local minimizers of $\Gamma$. To this end, let $f:\mathbb {R}\to \mathbb {R}$ be a continuous function satisfying

1. (f ₁). $f$ is locally Hölder continuous on $\mathbb {R}\backslash \{0\}$, $f(s)\equiv 0$ for $s\le 0$, and $f$ is strictly increasing in $(0,\,+\infty )$;

2. (f ₂). there exists some positive number $r_0\in (0,\,1)$ such that

   \[ \int_{0}^{s}f(t)\,{\rm d}t \le r_0 f(s)s, \quad \forall s\ge 0. \]

3. (f ₃). For all $\tau >0$,

   \[ \lim_{s\to +\infty}\left(f(s)\,e^{-\tau s}\right)=0. \]

Note that many nonlinearities that frequently appear in nonlinear elliptic equations satisfy ($f_1$)–($f_3$), for instance $f(s)=s_+^p$ with $p\in (0,\,+\infty )$. Our second result is as follows.

Theorem 1.3 Suppose that $k(\cdot )$ satisfies $(\mathcal {K})$ and $f$ satisfies $(\text {f}_1)$–$(\text {f}_3)$. Let $x_0$ be a strict local minimizer of $\Gamma$. Then there exists $\lambda _0>0$, such that for any $\lambda \in (\lambda _0,\,+\infty )$, (1.1) has a weak solution $(u^\lambda,\, c^\lambda )$ which satisfies properties as follows$:$

1. (1) the diameter of the plasma set $\{x\in D\mid u^\lambda (x)>0\}$ is of the order $O(\lambda ^{-({1}/{2})})$ as $\lambda \to +\infty$.

2. (2) For any $x\in \{x\in D\mid u^\lambda (x)>0\}$, $x$ tends to $x_0$ as $\lambda \to +\infty$.

3. (3) For $\lambda$ sufficiently large, $\{x\in D\mid u^\lambda (x)=0\}$ is a $C^1$ curve and converges to a circle as $\lambda \to +\infty$.

4. (4) There holds

   (1.10)\begin{equation} c^\lambda={-}\frac{I}{4\pi}\ln\lambda- \frac{I}{4\pi}\ln k(x_0)+Ih(x_0,x_0)-C_*+o(1). \end{equation}
   Here $C_*=({1}/{2\pi })\int _{\mathbb {R}^2}\ln ({1}/{|x^*-y'|})f(U)(y')\,{\rm d}y'$, where $U$ is the unique radial function satisfying
   (1.11)\begin{equation} \begin{cases} -\Delta U(x)= f(U)(x),\quad x\in \mathbb{R}^2,\\ \displaystyle\int_{\mathbb{R}^2}f(U)(x)\,{\rm d}x=I, \end{cases} \end{equation}
   and $x^*$ is any point of $\partial \{x\in \mathbb {R}^2\mid f(U)(x)=0\}$.

Remark 1.4 We give some examples to show the existence of $U$ in (1.11). If we choose $f(t)=t_+$, which corresponds to the classical confined plasma problem, then the unique radial $C^1$ solution of (1.11) has the explicit profile

\[ U(x)= \begin{cases} c\varphi_1(x), & |x|\leq s;\\ cs\varphi'_1(x)\ln\frac{|x|}{s}, & |x|\geq s, \end{cases} \]

where $s$ is a constant such that $1$ is the first eigenvalue of $-\Delta$ in $B_s(0)$ with Dirichlet boundary condition, $\varphi _1> 0$ is the first eigenfunction of $-\Delta$ in $B_s(0)$ with $\varphi _1(0) = 1$, $\varphi '_1$ is the derivative of $\varphi _1$, and $c$ is a constant such that $c\int _{B_s(0) }\varphi _1\,{\rm d}x=I.$ Note that $U\in C^{2,\alpha }(\mathbb {R}^2)$ for any $\alpha \in (0,\,1).$ For more results, see [9, 24] for example.

The idea of the proof of theorems 1.1 and 1.3 is the Arnold's variational principle, that is, considering maximization of some functional for the vorticity and analysing asymptotic behaviour of solutions. To this end, we introduce the definition of ‘vorticity’ $w=-\Delta u$, which originally comes from the study of the incompressible Euler equation, see [37]. Then we deduce the vorticity formulation (2.3) of (1.1) and give an equivalent description of main results, i.e., theorems 2.1 and 2.2. Indeed, we can generalize (2.3) to equation (2.5), which corresponds to solutions concentrating near several points. It suffices to prove the existence of solutions to (2.5) concentrating near strict local minimizers of an auxiliary function $\Gamma _l$, i.e., theorems 2.3 and 2.4, which is a generalized version of theorems 2.1 and 2.2. For the proof of theorem 2.3, the key is to regard the non-autonomous term $k$ as a $measure$. Note that because of the presence of the measure $k(x)\,{\rm d}x$ in the energy functional and the admissible class, the classical computation of vorticity method fails and we must give new estimates of maximizers, such as the energy $E$, the Lagrange multiplier $\mu ^\varepsilon$, the diameter and limiting location of the plasma set of $\omega ^\varepsilon$. For the proof of theorem 2.4, the differences from theorem 2.3 are as follows. First, to show the existence and profile of maximizers, we introduce another parameter $T$. Then we need to compute the upper bound of the stream function $\Psi _i^{\varepsilon,T}$ to eliminate the patch part and show that maximizers are solutions of (2.5), see lemma 4.5. Second, in order to get asymptotic behaviour of solutions, the limits of $\omega _i^{\varepsilon,T}$ and $\Psi _i^{\varepsilon,T}$ need to be estimated accurately.

Remark 1.5 We give some comments about the relation between our results and results in [12–14, 19, 24, 25, 28]. Note that in [12–14, 28], the term $k$ has to be a constant and the nonlinearity is $f(t)=t^p_+$ for $p\geq 0$. The key of proof is the use of the Lyapunov–Schmidt reduction method and the non-degeneracy of solutions to

(1.12)\begin{equation} -\Delta u=u^p_+,\quad \text{in}\ \mathbb{R}^2 \end{equation}

for $p\geq 0.$ Especially in [14], the nonlinearity $f$ is a Heaviside function and not differentiable and thus the proof requires very delicate estimates. Compared to these results, in this paper we can construct solutions of (1.1) with $k$ not a constant and general nonlinearity $f$. The key of proof is to use the expansion of Green's function $G(x,\,y)$ to prove the radius of the plasma set, the concentration location and the order of energy of the solution as $\lambda \to +\infty$. Indeed, the advantage of using Arnold's variational principle is that we do not need the non-degeneracy of solutions to (1.12) with $f(t)=t^p_+$ replaced by general $f$, which is also not known for general $f$. The argument adopted here is not affected by this issue and all we need is that the nonlinear term $f$ satisfies some growth conditions. This is why our result holds for general $f$. When choosing $f(t)=e^t$, del Pino et al. [19] constructed concentrated solutions of the equation

\[ -\Delta u=\varepsilon^2 K(x)e^u, \ \text{in}\ D;\quad u=0,\ \text{on}\ \partial D \]

such that the energy concentrates near small neighbourhoods of points $x_{1,\varepsilon },\,\ldots,\, x_{n,\varepsilon }$ as $\varepsilon \to 0$. These points tend to a critical point of the function defined by (1.7). When $n=1$, the function is $\Gamma$. Thus to some extent, the limiting behaviour of solutions in theorems 1.1 and 1.3 coincides with that in [19], and the only difference between theorems 1.1, 1.3 and results in [19] is the choice of the nonlinearity $f$. Finally, [24, 25] considered solutions of equation (1.5) in the case of $N\geq 3$ and $k$ not a constant. As $\varepsilon \to 0$, the plasma region of solutions to (1.5) will shrink to maximiers of $k$, rather than critical points of $\Gamma$. Note that the total vorticity vanishes as $\varepsilon \to 0$, that is, $\int _{\partial D} ({\partial u_\varepsilon }/{\partial \nu })\,{\rm d}s\to 0\ \text {as}\ \varepsilon \to 0.$ In contrast to these results, our result holds for $N=2$, and the total vorticity tends to a non-zero prescribed constant $I$ as $\varepsilon \to 0$. That is one of main differences between our results and results in [24, 25]. For more related works, see [4–6, 8, 20, 21, 39] for instance.

This paper is organized as follows. In § 2, we deduce the vorticity formulation of (1.1) and generalize main results to theorems 2.3 and 2.4, respectively. In § 3, we prove theorem 2.3 by solving a maximization problem of an energy functional for vorticity over admissible sets and giving asymptotic estimates of maximizers for $\varepsilon$ sufficiently small. The proof of theorem 2.4 will be shown in § 4.

2. Equivalent problem of (1.1)

We first reduce (1.1) to a dual problem for the vorticity. Let us define the vorticity $w=-\Delta u$. Since $u$ is a constant on $\partial D$, we have

\[ u(x)=\mathcal{G}w(x)-\mu=\int_{D}G(x,y)w(y)\,{\rm d}y-\mu,\quad x\in D \]

for some constant $\mu$, where $G(x,\,y)$ is the Green's function of $-\Delta$ in $D$ with zero Dirichlet condition. Taking this into (1.1) we have

(2.1)\begin{equation} w=\lambda k(x)f(\mathcal{G}w-\mu)\quad x\in D. \end{equation}

Using Green's formula, the third equation of (1.1) becomes

(2.2)\begin{equation} I={-}\int_{D}\frac{\partial u}{\partial \nu}\,{\rm d}s=\int_{D}w\,{\rm d}x. \end{equation}

Let us define $\varepsilon =\lambda ^{-({1}/{2})}$ and $\omega =w/k(x)$. Taking $\omega$ into (2.1) and (2.2), we get equations for $\omega$

(2.3)\begin{equation} \begin{cases} \displaystyle\omega=\frac{1}{\varepsilon^2} f(\mathcal{G}(k(x) \omega)-\mu), & x\in D,\\ \displaystyle\int_{D}\omega(x) k(x)\,{\rm d}x=I. \end{cases} \end{equation}

Note that it is equivalent to solve solution pairs $(u^\lambda,\, c^\lambda )$ of (1.1) and solution pairs $(\omega ^\varepsilon,\, \mu ^\varepsilon )$ of (2.3). Indeed, for a solution pair $(\omega ^\varepsilon,\, \mu ^\varepsilon )$ of (2.3), one can recover solutions of (1.1) by letting $u^\lambda =G(k(\cdot )\omega ^\varepsilon )-\mu ^\varepsilon$ and $c^\lambda =-\mu ^\varepsilon$.

For equation (2.3), we get the following equivalent description of theorems 1.1 and 1.3. Note that $\{x\in D\mid u^\lambda (x)>0\}=supp(\omega ^\varepsilon )$.

Theorem 2.1 Suppose that $k(\cdot )$ satisfies $(\mathcal {K})$ and $f(t)=\chi _{\{t>0\}}$. Let $x_0$ be a strict local minimizer of $\Gamma$. Then there exists $\varepsilon _0>0$, such that for any $\varepsilon \in (0,\,\varepsilon _0)$, (2.3) has a weak solution pair $(\omega ^\varepsilon,\, \mu ^\varepsilon )$ with the following properties:

1. (1) $diam(supp(\omega ^\varepsilon ))=O(\varepsilon )$.

2. (2) For any $x\in supp(\omega ^\varepsilon )$, $x$ tends to $x_0$ as $\varepsilon \to 0^+$.

3. (3) For $\varepsilon$ sufficiently small, $\partial (supp(\omega ^\varepsilon ))$ is a $C^1$ curve and converges to a circle as $\varepsilon \to 0^+$.

4. (4) There holds

   \[ \mu^{\varepsilon}= \frac{I}{2\pi}\ln\frac{1}{\varepsilon}+ \frac{k(x_0)}{2\pi}\int_{B_{\sqrt{{I}/{\pi k(x_0)}}}(0)}\ln\frac{1}{|x^*-y|}\,{\rm d}y-Ih(x_0,x_0)+o(1), \]
   where $x^*$ is any point of $\partial B_{\sqrt {{I}/{\pi k(x_0)}}}(0)$.

Theorem 2.2 Suppose that $k(\cdot )$ satisfies $(\mathcal {K})$ and $f$ satisfies $(\text {f}_1)$–$(\text {f}_3)$. Let $x_0$ be a strict local minimizer of $\Gamma$. Then there exists $\varepsilon _0>0$, such that for any $\varepsilon \in (0,\,\varepsilon _0)$, (2.3) has a weak solution pair $(\omega ^\varepsilon,\, \sigma ^\varepsilon )$ with the following properties:

1. (1) $diam(supp(\omega ^\varepsilon ))=O(\varepsilon )$.

2. (2) For any $x\in supp(\omega ^\varepsilon )$, $x$ tends to $x_0$ as $\varepsilon \to 0^+$.

3. (3) For $\varepsilon$ sufficiently small, $\partial (supp(\omega ^\varepsilon ))$ is a $C^1$ curve and converges to a circle as $\varepsilon \to 0^+$.

4. (4) There holds

   \[ \sigma^{\varepsilon}= \frac{I}{2\pi}\ln\frac{1}{\varepsilon}+ \frac{I}{4\pi}\ln k(x_0)-Ih(x_0,x_0)+C_*+o(1), \]
   where $C_*=({1}/{2\pi })\int _{\mathbb {R}^2}\ln ({1}/{|x^*-y'|})f(U)(y')\,{\rm d}y'$, $U$ is the unique radial function satisfying (1.11) and $x^*$ is any point of $\partial \{x\in \mathbb {R}^2\mid f(U)(x)=0\}$.

Indeed, to prove theorems 2.1 and 2.2, one can directly consider solutions of (1.1) concentrating near several distinct points. Let $l$ be an integer and $d_i\in \mathbb {R}/\{0\} (i=1,\,\ldots,\,l)$ be $l$ constants. Let us introduce an auxiliary function $\Gamma _l$ which generalizes the function $\Gamma$ defined by (1.8). Define

(2.4)\begin{equation} \Gamma_l(x_1,\ldots,x_{l})=\mathcal{H}_{l}(x_1,\ldots,x_{l})-\frac{1}{8\pi}\sum_{i=1}^ld_i^2\ln k(x_i), \end{equation}

where $\mathcal {H}_{l}$ is defined by (1.4). Notice that if $l=1$, then the auxiliary function is $\Gamma (x)$.

Let $(x_{0,1},\,\ldots,\,x_{0,l})$ be a strict local minimizer of $\Gamma _l$, that is, $(x_{0,1},\,\ldots,\,x_{0,l})$ is the unique minimizer of $\Gamma _l$ over $\overline {B_1}\times \cdots \times \overline {B_l}$. Here $B_i:=B_{\delta }(x_{0,i})$ for some $\delta >0$ sufficiently small such that $\overline {B_i}\subset D$ and $\overline {B_i}\cap \overline {B_j}=\varnothing$ for $i\neq j$. Consider solution pairs $(\omega ^\varepsilon,\, \mu ^\varepsilon _i) (i=1,\,\ldots,\,l)$ of the following equations

(2.5)\begin{equation} \begin{cases} \displaystyle\omega=\sum_{i=1}^l\frac{sgn(d_i)}{\varepsilon^2}f_i(\mathcal{G}(k \omega)-\mu_i)\chi_{B_i}, & x\in D,\\ \displaystyle\int_{B_i}\omega k(x)\,{\rm d}x=d_i, \end{cases} \end{equation}

where $sgn(d_i)=1$ if $d_i>0$ and $sgn(d_i)=-1$ if $d_i<0$. $f_i$ are $l$ given functions and $\mu _i$ are unknown constants. The following result shows that for any strict local minimizer $(x_{0,1},\,\ldots,\,x_{0,l})$ of $\Gamma _l$, there exist solutions of (2.5) concentrating near $l$ distinct points $x_{0,i}$.

Theorem 2.3 Suppose that $k$ satisfies $(\mathcal {K})$ and $f_i(t)=\chi _{\{t>0\}}$ ($i=1,\,\ldots,\,l$). Then for any strict local minimizer $(x_{0,1},\,\ldots,\,x_{0,l})$ of $\Gamma _l$, there exists $\varepsilon _0>0$, such that for any $\varepsilon \in (0,\,\varepsilon _0)$, (2.5) has a weak solution pair $(\omega ^\varepsilon,\, \mu ^\varepsilon _i)$ with the following properties:

1. (1) Define $\omega _i=\omega \chi _{ B_i}$. Then $diam(supp(\omega _i^\varepsilon ))=O(\varepsilon )$.

2. (2) The support of $\omega _i^\varepsilon$ tends to $x_{0,i}$ as $\varepsilon \to 0^+$, that is,

   \[ \lim_{\varepsilon \to 0^+}\sup_{x\in supp(\omega_i^\varepsilon)}|x- x_{0,i}|=0. \]

3. (3) For $\varepsilon$ sufficiently small, $\partial (supp(\omega _i^\varepsilon ))$ is a $C^1$ curve and converges to a circle as $\varepsilon \to 0^+$.

4. (4) There holds

   \begin{align*} \mu_i^{\varepsilon}& = \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+ \frac{k(x_{0,i})}{2\pi}\int_{B_{\sqrt{{|d_i|}/{\pi k(x_{0,i})}}}(0)}\ln\frac{1}{|x^*-y|}\,{\rm d}y-|d_i|h(x_{0,i},x_{0,i})\\ & \quad+sgn(d_i)\sum_{j=1,j\neq i}^ld_jG(x_{0,i},x_{0,j})+o(1), \end{align*}
   where $x^*$ is any point of $\partial B_{B_{\sqrt {{|d_i|}/{\pi k(x_{0,i})}}}(0)}$.

Theorem 2.4 Suppose that $k(\cdot )$ satisfies $(\mathcal {K})$ and $f_i$ satisfies $(\text {f}_1)-(\text {f}_3)$ ($i=1,\,\ldots,\,l$). Then for any strict local minimizer $(x_{0,1},\,\ldots,\,x_{0,l})$ of $\Gamma _l$, there exists $\varepsilon _0>0$, such that for any $\varepsilon \in (0,\,\varepsilon _0)$, (2.5) has a weak solution pair $(\omega ^\varepsilon,\, \sigma ^\varepsilon _i)$ with the same properties as those in theorem 2.3. Moreover, $\sigma ^\varepsilon _i$ has the following estimates

\begin{align*} \sigma_i^{\varepsilon}& = \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+ \frac{|d_i|}{4\pi}\ln k(x_{0,i})-|d_i|h(x_{0,i},x_{0,i})\\ & \quad +sgn(d_i)\sum_{j=1,j\neq i}^ld_jG(x_{0,i},x_{0,j})+C_i+o(1). \end{align*}

Here $C_i=({1}/{2\pi })\int _{\mathbb {R}^2}\ln ({1}/{|x^*-y'|})f_i\circ U_i(y')\,{\rm d}y'$, where $U_i$ is the unique radial function satisfying

\[ \begin{cases} -\Delta U_i(x)= f_i(U_i)(x),\quad x\in \mathbb{R}^2,\\ \int_{\mathbb{R}^2}f_i(U_i)(x)\,{\rm d}x=|d_i|, \end{cases} \]

and $x^*$ is any point of $\partial \{x\in \mathbb {R}^2\mid f_i\circ U_i(x)=0\}$.

Note that when $l=1$, results in theorems 2.3 and 2.4 reduce to those in theorems 2.1 and 2.2. We will give the proof of theorems 2.3 and 2.4 directly in the following sections.

3. Proof of theorem 2.3

3.1. Variational problem

In this section, we give the proof of theorem 2.3. We define the energy functional

(3.1)\begin{equation} E(\omega)=\frac{1}{2}\int_{D}\int_{D}G(x,y)\omega(x)\omega(y)\,{\rm d}\theta(x)\,{\rm d}\theta(y), \end{equation}

where $\,{\rm d}\theta (x)=k(x)\,{\rm d}x$ is a $measure$ deduced by the non-autonomous term $k$. By the assumption $(\mathcal {K})$, we know that $\,{\rm d}\theta (x)$ is equivalent to the two-dimensional Lebesgue measure $\,{\rm d}x$.

Define a constraint set

(3.2)\begin{align} \mathcal{M}_{\varepsilon}(D)& =\left\{\omega=\sum_{i=1}^{l}\omega_i\in L^\infty(D)\mid\ \omega_i=\omega \chi_{B_i},\right. \notag\\ & \qquad \left. 0\le sgn(d_i)\omega_i \le \frac{1}{ \varepsilon^2},\ \int_{B_i}\omega_i\,{\rm d}\theta(x) =d_i\vphantom{\sum_{i=1}^{l}}\right\}. \end{align}

The difference between $\mathcal {M}_{\varepsilon }(D)$ and the classical results is that we impose the $L^1(B_i,\, \,{\rm d}\theta (x))$ norm of $\omega _i$ to be $d_i$, rather than the $L^1(B_i,\, \,{\rm d}x)$ norm, which may cause essential difficulty in proving asymptotic behaviour of solutions.

Consider the maximization problem

\[ (\mathcal{P}')\quad \sup_{\omega\in \mathcal{M}_\varepsilon(D)}E(\omega). \]

To begin with, we show the existence and profile of maximizers of $E$ over $\mathcal {M}_{\varepsilon }(D)$.

Proposition 3.1 There exists $\omega ^{\varepsilon }\in \mathcal {M}_{\varepsilon }(D)$, such that $E(\omega ^{\varepsilon })=\sup _{\tilde {\omega }\in \mathcal {M}_{\varepsilon }(D)}E(\tilde {\omega })$.

Proof. Since $G(\cdot,\,\cdot )\in L^1(D\times D)$, we know that $E$ is bounded from above on the set $\mathcal {M}_{\varepsilon }(D)$. Now we choose a maximization sequence $\{\omega ^n\}\subset \mathcal {M}_{\varepsilon }(D)$ of $E$, that is,

\[ \lim_{k\to +\infty}E(\omega^n)=\sup_{\omega \in \mathcal{M}_{\varepsilon}(D)}E(\omega). \]

By direct computations we can prove that $\mathcal {M}_{\varepsilon }$ is a sequentially compact subset of $L^2(D)$ in the weak topology. So we may assume that, up to a subsequence, $\omega ^n\to \omega ^{\varepsilon }$ weakly in $L^2(D)$ as $n\to +\infty$ for some $\omega ^{\varepsilon }\in \mathcal {M}_{\varepsilon }$. So

\[ k(x)\omega^n\to k(x)\omega^{\varepsilon } \quad \text{in}\ L^2(D) \text{weak topology}. \]

By elliptic regularity theory,

\[ \mathcal{G}(k\omega^{n})=\int_{D}G(x,y)\omega^{n}\,{\rm d}\theta(x)\to \int_{D}G(x,y)\omega^{\varepsilon}\,{\rm d}\theta(x)=\mathcal{G}(k \omega^{\varepsilon })\quad \text{in}\ W^{1,p}(D) \]

for any $p>1$, from which we deduce that

\[ \lim_{n\to+\infty}E(\omega^n) =E(\omega^\varepsilon). \]

So $\omega ^{\varepsilon }$ is a maximizer of $E$ over $\mathcal {M}_{\varepsilon }$.

We define $\omega _i^{\varepsilon }=\omega ^{\varepsilon } \chi _{B_i}$ to be each piece of the maximizer. Then using classical idea in [37] we can get that the maximizers has the form of (2.5).

Proposition 3.2 Let $\omega ^{\varepsilon }$ be a maximizer defined as in lemma 3.1. Then

(3.3)\begin{equation} \omega^\varepsilon=\sum_{i=1}^l\omega_i^\varepsilon =\sum_{i=1}^l\frac{sgn(d_i)}{ \varepsilon^2} \chi_{\left\{\psi_i^\varepsilon>0\right\}\cap B_i}, \end{equation}

where $\psi _i^\varepsilon :=sgn(d_i)\mathcal {G}(k \omega ^\varepsilon )-\mu ^\varepsilon _i$ and $\mu ^\varepsilon _i$ are constants dependent on $\varepsilon$ for $i=1,\,\ldots,\,l$. Moreover, for $\varepsilon$ sufficiently small there holds

(3.4)\begin{equation} \mu_i^{\varepsilon}\geq C_0, \end{equation}

where $C_0$ is a negative constant independent of $\varepsilon$.

Proof. Let $\omega ^{\varepsilon }$ be a maximizer. For any ${\omega }\in \mathcal {M}_{\varepsilon }$, we set

\[ \omega_{(s)}=\omega^{\varepsilon }+s({\omega}-\omega^{\varepsilon }),\quad \text{for} \ s\in[0,1]. \]

Since $\mathcal {M}_{\varepsilon }$ is a convex set, $\omega _{(s)}\in \mathcal {M}_{\varepsilon }$ for any $s\in [0,\,1]$. So $E(\omega _{(s)})\leq E(\omega ^\varepsilon )$, which implies that

\[ 0 \ge \left.\frac{{\rm d} E(\omega_{(s)})}{{\rm d}s}\right|_{s=0^+} =\int_{D}({\omega}-\omega^{\varepsilon}) \mathcal{G}(k \omega^{\varepsilon}) \,{\rm d}\theta(x), \]

that is,

\[ \int_{D}\omega \mathcal{G}(k \omega^{\varepsilon}) \,{\rm d}\theta(x) \leq \int_{D}\omega^{\varepsilon } \mathcal{G}(k \omega^{\varepsilon}) \,{\rm d}\theta(x) \]

for any ${\omega }\in \mathcal {M}_{\varepsilon }(D).$ By the definition of $\mathcal {M}_{\varepsilon }(D)$ and the bathtub principle (see [26]), we get for any $i=1,\,\ldots,\,l$

(3.5)\begin{equation} \begin{aligned} sgn(d_i) \mathcal{G}(k \omega^{\varepsilon}) & \ge \mu_i^{\varepsilon} \quad \mbox{on} \left\{sgn(d_i)\omega_i^{\varepsilon}=\frac{1}{\varepsilon^2}\right\}\cap B_i, \\ sgn(d_i) \mathcal{G}(k \omega^{\varepsilon}) & = \mu_i^{\varepsilon} \quad \mbox{on}\ \left\{0< sgn(d_i)\omega_i^{\varepsilon}<\frac{1}{\varepsilon^2}\right\}\cap B_i, \\ sgn(d_i) \mathcal{G}(k \omega^{\varepsilon}) & \le \mu_i^{\varepsilon} \quad \mbox{on} \{sgn(d_i)\omega_i^{\varepsilon}=0\}\cap B_i, \end{aligned} \end{equation}

where $\mu _i^{\varepsilon }$ is a constant satisfying

(3.6)\begin{equation} \mu_i^{\varepsilon}=\inf\{s\in\mathbb{R}\mid|\{x\in B_i\mid sgn(d_i) \mathcal{G}\left(k \omega^{\varepsilon} \right) >s\}|_{\theta}\le |d_i|\varepsilon^2\}. \end{equation}

Notice that $|\{x\in B_i\mid sgn(d_i) \mathcal {G}(k \omega ^{\varepsilon } ) >s\}|_{\theta }$ means that the $\,{\rm d}\theta (x)-$measure of the set $\{x\in B_i\mid sgn(d_i) \mathcal {G}(k \omega ^{\varepsilon } ) >s\}$. Thus using (3.5), we get that

\[ \omega_i^{\varepsilon}= sgn(d_i)\frac{1}{ \varepsilon^2}\chi_{\{ \psi_i^{\varepsilon}>0 \}\cap B_i }, \]

where $\psi _i^{\varepsilon }=sgn(d_i)\mathcal {G}(k \omega ^{\varepsilon } )-\mu ^{\varepsilon }_i$. So (3.3) is proved.

It remains to prove that $\mu _i^{\varepsilon }\geq C_0$ for some $C_0$. In fact for any $x\in B_i$,

(3.7)\begin{align} sgn(d_i)\mathcal{G}\left( k \omega^{\varepsilon } \right)(x) & =\mathcal{G}(k |\omega_i^{\varepsilon }|)(x)+sgn(d_i)\sum_{j\neq i}\mathcal{G}(k \omega_j^{\varepsilon })(x) \nonumber\\ & \geq{-}\sum_{j\neq i}k_1|d_j|\max_{x\in B_i, y\in B_j}|G(x,y)|. \end{align}

Combining (3.6) with (3.7), we can get (3.4) by letting $C_0=-\sum _{j\neq i}k_1|d_j| \max _{x\in B_i, y\in B_j}|G(x,\,y)|$.

Remark 3.3 Indeed, one can repeat the proof in [37] to prove proposition 3.2. Without loss of generality, we assume $d_i>0$. For any $z_1,\,z_2 \in L^{\infty }(D)$ satisfying

(3.8)\begin{equation} \begin{cases} supp(z_1),supp(z_2)\subset B_i,\\ \displaystyle z_1,z_2\geq 0,\,\,\text{a.e. in } D, \int_{D}z_1(x)\,{\rm d}\theta(x)=\int_{D}z_2(x)\,{\rm d}\theta(x),\\ \displaystyle z_1=0 \quad \text{in }\ D\setminus\left\{x\in D\mid\omega^{\varepsilon}(x)\leq \frac{1}{ \varepsilon^2}-a\right\}, \\ \displaystyle z_2=0 \quad \text{in } D\setminus\left\{x\in D\mid\omega^{\varepsilon}(x)\geq a\right\}, \end{cases} \end{equation}

where $a>0$ is sufficiently small, we define a family of functions $\omega _s=\omega ^{\varepsilon }+ s(z_1-z_2)$, $s>0$. Then one can prove that $\omega _s\in M_{\varepsilon }(D)$ for $s>0$ sufficiently small. So ${{\rm d}E(\omega _s)}/{{\rm d}s}|_{s=0^+}\leq 0$, which implies that

\[ \int_{D}\mathcal{G}(k \omega^{\varepsilon})\left(x\right)z_1(x)\,{\rm d}\theta(x)\leq \int_{D}\mathcal{G}(k \omega^{\varepsilon})\left(x\right)z_2(x)\,{\rm d}\theta(x). \]

From this we get

\[ \sup_{\left\{x\in D|\omega^{\varepsilon}(x)< {1}/{\varepsilon^2}\right\}\cap B_i} \mathcal{G}(k \omega^{\varepsilon})\left(x\right)= \inf_{\left\{x\in D|\omega^{\varepsilon}(x)>0\right\}\cap B_i} \mathcal{G}(k \omega^{\varepsilon})\left(x\right). \]

Define $\mu _i^{\varepsilon }:=\inf _{\left \{x\in D|\omega ^{\varepsilon }(x)>0\right \}\cap B_i} \mathcal {G}(k \omega ^{\varepsilon })(x),$ it is not hard to prove that

\[ \begin{cases} \displaystyle\omega^{\varepsilon}=\frac{1}{ \varepsilon^2} \quad \text{a.e. in }\ B_i\cap \left\{x\in D\mid \mathcal{G}(k \omega^{\varepsilon})>\mu_i^{\varepsilon}\right\}, \\ \omega^{\varepsilon}=0 \quad \text{a.e. in }\ B_i\cap \left\{x\in D\mid \mathcal{G}(k \omega^{\varepsilon})<\mu_i^{\varepsilon}\right\}. \end{cases} \]

On $\left \{x \in D\mid \mathcal {G}(k \omega ^{\varepsilon })=\mu _i^{\varepsilon }\right \}$, by properties of Sobolev space, we have $k \omega ^{\varepsilon }=0$. So $\omega ^{\varepsilon }=0$ a.e. in $\left \{x \in D\mid \mathcal {G}(k \omega ^{\varepsilon })=\mu _i^{\varepsilon }\right \}$. Thus

\[ \omega_i^{\varepsilon}=\frac{1}{ \varepsilon^2}\chi_{\{x\in D\mid \mathcal{G}(k \omega^{\varepsilon})>\mu_i^{\varepsilon} \}\cap B_i }. \]

So using this method, we can also get the same results as that in proposition 3.2.

3.2. Asymptotic analysis

In the following, we give asymptotic estimates of $\omega ^\varepsilon _i$. We first give lower bound of the energy $E_i(\omega ^{\varepsilon })$ and the Lagrange multiplier $\mu _i^\varepsilon$. Note that since the measure in (3.1) is $k(x)\,{\rm d}x$, we need to choose test functions properly. Then using the properties of function $\ln x$ and the theory of rearrangement function, we get that the diameter of the plasma set of $\omega ^{\varepsilon }$ is the order of $\varepsilon$ and the limiting location is a minimizer of $\Gamma _l$.

To simplify the proof , we define the energy functional associated with $\omega ^\varepsilon _i$

(3.9)\begin{equation} E_i(\omega):= \frac{1}{2}\int_{D}G(x,y)\omega_i(x)\omega_i(y) \,{\rm d}\theta(x)\,{\rm d}\theta(y)\quad \text{for}\ i=1,\ldots,l. \end{equation}

Direct computation shows that

(3.10)\begin{equation} E(\omega)=\sum_{k=1}^lE_k(\omega)+O(1)=E_i(\omega)+E\left(\sum_{j\neq i}\omega_j\right)+O(1)\quad \forall \omega\in M_{\varepsilon}(D). \end{equation}

We first give a rough lower bound of $E_i(\omega ^{\varepsilon })$.

Lemma 3.4 Let $\omega ^{\varepsilon }$ be a maximizer. Then for $i=1,\,\ldots,\,l$

\[ E_i(\omega^{\varepsilon})\geq \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+O(1). \]

Proof. We choose

\[ \tilde{\omega}^{\varepsilon}=\frac{sgn(d_i)}{\varepsilon^2}\chi_{B_{t(\varepsilon)\varepsilon}(x_{0,i})}+\sum_{k=1,k\neq i}^{l}\omega_k^{\varepsilon}=\sum_{k=1}^{l}\tilde{\omega}_k^{\varepsilon}, \]

where $t(\varepsilon )\in (({1}/{2})\sqrt {{|d_i|}/{\pi k(x_{0,i})}},\, 2\sqrt {{|d_i|}/{\pi k(x_{0,i})}})$ satisfies $\int _{B_i}\tilde {\omega }_i^{\varepsilon }\,{\rm d}\theta (x)=d_i$. Direct calculations show that $t(\varepsilon )$ exists for $\varepsilon$ sufficiently small and $\lim _{\varepsilon \to 0^+}t(\varepsilon )=\sqrt {{|d_i|}/{\pi k(x_{0,i})}}$. Then $\tilde {\omega }^{\varepsilon }\in \mathcal {M}_{\varepsilon }(D)$ and $E(\omega ^\varepsilon )\geq E(\tilde {\omega }^{\varepsilon })$, which implies that

(3.11)\begin{align} E(\omega^\varepsilon)& \geq{-}\frac{1}{4\pi}\int_{D}\int_{D}\ln|{x}-{y}|\tilde{\omega}_i^\varepsilon({x})\tilde{\omega}_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y)\notag\\ & \quad - \frac{1}{2}\int_{D}\int_{D}h({x},{y})\tilde{\omega}_i^\varepsilon({x})\tilde{\omega}_i^\varepsilon({y})\theta(x)\,{\rm d}\theta(y) \nonumber\\ & \quad +E\left(\sum_{j\neq i}\omega_j^{\varepsilon}\right)+O(1). \end{align}

Since the diameter of $supp(\tilde {\omega }_i^\varepsilon )$ is $\sqrt {({|d_i|}/{\pi k(x_{0,i})})}\varepsilon +o(\varepsilon )$, we have

(3.12)\begin{equation} -\frac{1}{4\pi}\int_{D}\int_{D}\ln|{x}-{y}|\tilde{\omega}_i^\varepsilon({x})\tilde{\omega}_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y) \geq \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+O(1). \end{equation}

By the choice of $B_i$ we obtain

(3.13)\begin{equation} \left|\frac{1}{2}\int_{D}\int_{D}h({x},{y})\tilde{\omega}_i^\varepsilon({x})\tilde{\omega}_i^\varepsilon({y})\theta(x)\,{\rm d}\theta(y)\right| = O(1), \end{equation}

Taking (3.12), (3.13) into (3.11) and using (3.10), we get the desired result.

Then, one can get the lower bound of Lagrange multipliers $\mu _i^\varepsilon$.

Lemma 3.5 Let $\omega ^{\varepsilon }$ be a maximizer and $\mu _i^\varepsilon$ be the associated Lagrange multiplier. Then there holds

(3.14)\begin{equation} \mu_i^{\varepsilon}\geq \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+O(1). \end{equation}

Proof. Let us first prove that

(3.15)\begin{equation} sgn(d_i)\int_{D}(sgn(d_i)\mathcal{G}(k\omega^\varepsilon)-\mu_i^{\varepsilon})\omega_i^\varepsilon \,{\rm d}\theta(x)=O(1). \end{equation}

Using the definition of $B_i$ we get

(3.16)\begin{align} & sgn(d_i)\int_{D}(sgn(d_i)\mathcal{G}(k\omega^\varepsilon)-\mu_i^{\varepsilon})\omega_i^\varepsilon \,{\rm d}\theta(x)\notag\\ & \quad =\int_{D}(sgn(d_i)\mathcal{G}(k\omega^\varepsilon)-\mu_i^{\varepsilon}+C_0)|\omega_i^\varepsilon| \,{\rm d}\theta(x)+|C_0d_i| \nonumber\\ & \quad \leq \int_{D}(sgn(d_i)\mathcal{G}(k\omega_i^\varepsilon)-\mu_i^{\varepsilon}+C_0)_+|\omega_i^\varepsilon| \,{\rm d}\theta(x)+ O(1). \end{align}

Define $P^\varepsilon _i=(sgn(d_i)\mathcal {G}(k\omega _i^\varepsilon )-\mu _i^{\varepsilon }+C_0)_+$ and $\bar {P}^\varepsilon _i=(sgn(d_i)\mathcal {G}(k\omega ^\varepsilon )-\mu _i^{\varepsilon }+C_0)_+$. So by (3.4), $P^\varepsilon _i,\,\bar {P}^\varepsilon _i\in H^1_0(D)$. On the one hand, we get

(3.17)\begin{equation} \int_{D}P^\varepsilon_i \omega_i^\varepsilon \,{\rm d}\theta(x)= \int_{D}|\nabla P^\varepsilon_i|^2 \,{\rm d}{x}. \end{equation}

On the other hand, by the choice of $C_0$ we have $supp(\bar {P}^\varepsilon _i)\cap B_i\subseteq supp(\omega _i^\varepsilon ),\,$ which implies that

\begin{align*} \int_{D}P^\varepsilon_i |\omega_i^\varepsilon| \,{\rm d}\theta(x) & \leq \int_{D}\bar{P}^\varepsilon_i |\omega_i^\varepsilon| \,{\rm d}\theta(x)+O(1)\\ & \leq\frac{1}{ \varepsilon^2}\int_{supp(\omega_i^\varepsilon)}\bar{P}^\varepsilon_i \,{\rm d}\theta(x)+O(1)\\ & \leq \frac{k_1}{ \varepsilon^2}|supp(\omega_i^\varepsilon)|^{{1}/{2}}\left(\int_{supp(\omega_i^\varepsilon)}(\bar{P}^\varepsilon_i)^2\,{\rm d}x\right)^{{1}/{2}}+O(1). \end{align*}

Using the continuous embedding $W^{1,1}(B_i)\subset L^2(B_i)$, we have

\begin{align*} & \int_{D}P^\varepsilon_i |\omega_i^\varepsilon| \,{\rm d}\theta(x)\\ & \quad \leq\frac{Ck_1}{ \varepsilon^2}|supp(\omega_i^\varepsilon)|^{{1}/{2}}\left(\int_{B_i}\bar{P}^\varepsilon_i\,{\rm d}{x}+\int_{B_i}|\nabla \bar{P}^\varepsilon_i| \,{\rm d}{x}\right)+O(1)\\ & \quad \leq \frac{Ck_1}{ \varepsilon^2}|supp(\omega_i^\varepsilon)|^{{1}/{2}}\left(\int_{supp(\omega_i^\varepsilon)}P^\varepsilon_i\,{\rm d}{x}+\int_{supp(\omega_i^\varepsilon)}|\nabla P^\varepsilon_i| \,{\rm d}{x}\right)+O(1)\\ & \quad \leq \frac{Ck_1}{k_0}|supp(\omega_i^\varepsilon)|^{{1}/{2}}\int_{D}P^\varepsilon_i\omega_i^\varepsilon \,{\rm d}\theta(x)\notag\\ & \qquad + \frac{Ck_1 }{ \varepsilon^2}|supp(\omega_i^\varepsilon)|^{{1}/{2}} \int_{supp(\omega_i^\varepsilon)}|\nabla P^\varepsilon_i| \,{\rm d}{x}+O(1). \end{align*}

So for $\varepsilon$ sufficiently small, we get $\int _{D}P^\varepsilon _i |\omega _i^\varepsilon | \,{\rm d}\theta (x)\leq ({Ck_1 }/{ \varepsilon ^2})|supp(\omega _i^\varepsilon )|^{{1}/{2}} \int _{supp(\omega _i^\varepsilon )} |\nabla P^\varepsilon _i| \,{\rm d}{x}+O(1).$ By Hölder's inequality,

(3.18)\begin{align} \int_{D}P^\varepsilon_i |\omega_i^\varepsilon| \,{\rm d}\theta(x)& \leq\frac{Ck_1 }{ \varepsilon^2}|supp(\omega_i^\varepsilon)| \left(\int_{supp(\omega_i^\varepsilon)}|\nabla P^\varepsilon_i|^2 \,{\rm d}{x}\right)^{{1}/{2}}+O(1) \nonumber\\ & \leq C k_1 \left(\int_{D}|\nabla P^\varepsilon_i|^2 \,{\rm d}{x}\right)^{{1}/{2}}+O(1). \end{align}

Combining (3.17), (3.18) and (3.16), we get (3.15).

Notice that

\[ 2E_i(\omega^\varepsilon)=sgn(d_i)\int_{D}(sgn(d_i)\mathcal{G}(k\omega^\varepsilon)-\mu_i^{\varepsilon})\omega_i^\varepsilon \,{\rm d}\theta(x)+|d_i|\mu^\varepsilon_i. \]

So using lemma 3.5, we get (3.14).

By lemma 3.5, one can prove that the diameter of $\omega _i^{\varepsilon }$ is $O(\varepsilon )$.

Lemma 3.6 Let $\omega ^{\varepsilon }$ be a maximizer. Then

\[ r_1\varepsilon\leq diam\left(supp( \omega_i^\varepsilon)\right)\le R_1\varepsilon, \]

where $r_1,\, R_1>0$ are constants independent of $\varepsilon$.

Proof. By the choice of $\mathcal {M}_\varepsilon (D)$, we know that $k_1|supp(\omega _i^\varepsilon )|\geq |d_i| \varepsilon ^2.$ This implies that

\[ diam\left(supp( \omega_i^\varepsilon)\right)\geq r_1\varepsilon \]

for some $r_1> 0$.

On the other hand, for any $x\in supp( \omega _i^\varepsilon )$, using (3.3) we have $\psi _i^{\varepsilon }(x)\geq 0$, which shows that

\begin{align*} \mu_i^{\varepsilon}& \leq{-}\frac{1}{2\pi}\int_{D}\ln|x-y ||\omega_i^{\varepsilon}|(y)\,{\rm d}\theta(y)-\int_{D}h(x,y)|\omega_i^\varepsilon(y)|\,{\rm d}\theta(y)\notag\\ & \quad +sgn(d_i)\sum_{j=1,j\neq i}^l\mathcal{G}(k\omega_j^{\varepsilon})(x)\\ & \leq{-}\frac{1}{2\pi}\int_{D}\ln|x-y||\omega_i^{\varepsilon}|(y)\,{\rm d}\theta(y)+O(1). \end{align*}

Thus by (3.14), we have

\[ \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}\leq{-}\frac{1}{2\pi}\int_{D}\ln|x-y||\omega_i^{\varepsilon}|(y)\,{\rm d}\theta(y)+O(1). \]

From the classical estimates in [37], we get $diam(supp( \omega _i^\varepsilon ))\le R_1\varepsilon$ for some $R_1>1$.

We now estimate the limiting location of $\omega _i^{\varepsilon }$ as $\varepsilon$ tends to 0. To begin with, we define the $\theta$-weighted mass centre of $\omega _i^\varepsilon$ as

\[ \bar{X}_i^\varepsilon:= \frac{1}{d_i}\int_{D}x \omega_{i}^\varepsilon(x) \,{\rm d}\theta(x)\ \ \text{ for}\ i=1,\ldots,l. \]

Then $\bar {X}_i^\varepsilon \in \overline {B_i}$. Since $\overline {B_i}$ is compact, we may choose a subsequence of $\{\bar {X}_i^{\varepsilon _n}\}_{n=1}^{ \infty }$(still denoted by $\bar {X}_i^\varepsilon$) satisfying

\[ \lim_{\varepsilon\to 0^+}\bar{X}_i^\varepsilon=\mathbf{x^*_i}\in \overline{B_i}. \]

Define the scaled function of $\omega _{i}^\varepsilon$

\[ \zeta^\varepsilon_i:= \frac{sgn(d_i)}{\varepsilon^2}\omega^\varepsilon_i(\varepsilon x+\bar{X}_i^\varepsilon) \ \ \ \ x\in D_\varepsilon, \]

where $D_\varepsilon =\{x\in \mathbb {R}^2\mid \varepsilon x+\bar {X}_i^\varepsilon \in D\}$. Then using the definition of $M_\varepsilon (D)$, we have $0\leq \zeta ^\varepsilon _i\leq 1.$ Moreover, by lemma 3.6 we get that the support set of $\zeta ^\varepsilon _i$ is contained in $B_{R_1}(0)$. Using $\int _{B_i}\omega _i^\varepsilon \,{\rm d}\theta (x)=d_i$, we get

(3.19)\begin{equation} \int_{B_{R_1}(0)}k(\varepsilon x+\bar{X}_i^\varepsilon)\zeta^\varepsilon_i(x)\,{\rm d}x=\int_{B_i}|\omega_i^\varepsilon|\,{\rm d}\theta(x)=|d_i|. \end{equation}

Since $\zeta ^\varepsilon _i$ is uniformly bounded in $L^p(B_{R_1}(0))$ for any $p\in [1,\,+\infty ]$, then still up to a subsequence, we may assume that $\zeta ^\varepsilon _i\to \zeta ^*_i$ in $L^p$ weak topology and $L^{\infty }$ weak star topology for some $\zeta ^*_i\in L^\infty (B_{R_1}(0))$ as $\varepsilon \to 0.$ We now calculate the necessary condition of $\mathbf {x^*_i}$ and the profile of $\zeta ^*_i$. To this end, we define a real-valued function

(3.20)\begin{equation} Q_i(t)=\frac{t^2}{4\pi}\iint\ln\frac{1}{|x-y|}\chi_{B_{\sqrt{{|d_i|}/{\pi t}}}(0)}(x)\chi_{B_{\sqrt{{|d_i|}/{\pi t}}}(0)}(y)\,{\rm d}x\,{\rm d}y. \end{equation}

Direct calculation shows that $Q_i(t)=({d_i^2}/{8\pi })\ln t+C^*$, where $C^*$ is a universal constant.

Proposition 3.7 There holds

(3.21)\begin{align} \mathcal{H}_l(\mathbf{x^*_1},\mathbf{x^*_2},\ldots,\mathbf{x^*_l})-\sum_{i=1}^lQ_i(k(\mathbf{x^*_i}))=\min_{x_i\in B_i}\mathcal{H}_l(x_1,x_2,\ldots,x_l)-\sum_{i=1}^lQ_i(k(x_i)). \end{align}

As a consequence, $\mathbf {x^*_i}=x_{0,i}$ for $i=1,\,\ldots,\,l.$ Moreover,

(3.22)\begin{equation} \zeta^\varepsilon_i\to\zeta^*_i=\chi_{B_{\sqrt{{|d_i|}/{\pi k(\mathbf{x^*_i})}}}(0)} \end{equation}

in $L^p$ topology for any $p>1$ as $\varepsilon \to 0.$

Proof. By proposition 3.2, we know that $|\omega ^\varepsilon _i|$ is a vortex patch with height ${1}/{\varepsilon ^2}$, so $\zeta ^\varepsilon _i$ is a vortex patch with height 1. So the limiting function $\zeta ^*_i$ is also a vortex patch with height 1, that is, $\zeta ^*_i=\chi _{U^*}$ for some set $U^*\subseteq B_{R_1}(0)$.

Since $k$ is a $C^0$ function and $\lim _{\varepsilon \to 0^+}\bar {X}_i^\varepsilon =\mathbf {x^*_i}$, we have

\[ k(\varepsilon x+\bar{X}_i^\varepsilon)\to k(\mathbf{x^*_i})\quad \text{uniformly in}\ B_{R_1}(0), \]

so by (3.19)

(3.23)\begin{equation} |d_i|=\lim_{\varepsilon \to 0^+}\int_{B_{R_1}(0)}k(\varepsilon x+\bar{X}_i^\varepsilon)\zeta^\varepsilon_i(x)\,{\rm d}x= k(\mathbf{x^*_i})|U^*|. \end{equation}

On the one hand, by the definition of $E(\omega ^\varepsilon )$ and lemma 3.6, we get

\begin{align*} E(\omega^\varepsilon)& ={-}\frac{1}{4\pi}\sum_{i=1}^l\iint\ln|{x}-{y}| \omega_i^\varepsilon({x}) \omega_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y)\\ & \quad - \frac{1}{2}\sum_{i=1}^l\iint h({x},{y})\omega_i^\varepsilon({x})\omega_i^\varepsilon({y})\theta(x)\,{\rm d}\theta(y) \\ & \quad +\frac{1}{2}\sum_{1\leq i\neq j\leq l}\iint G(x,y)\omega_i^\varepsilon({x})\omega_j^\varepsilon({y})\theta(x)\,{\rm d}\theta(y)\\ & = \sum_{i=1}^l\frac{1}{4\pi}\iint\ln\frac{1}{\varepsilon|{x}-{y}|} k(\varepsilon x+\bar{X}_i^\varepsilon)\zeta_i^\varepsilon({x}) k(\varepsilon y+\bar{X}_i^\varepsilon)\zeta_i^\varepsilon({y})\,{\rm d}x\,{\rm d}y\\ & \quad -\mathcal{H}_l(\mathbf{x^*_1},\mathbf{x^*_2},\ldots,\mathbf{x^*_l})+o(1)\\ & = \sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^l \frac{1}{4\pi}\int_{B_{R_1}(0)}\int_{B_{R_1}(0)}\\ & \quad \times \ln\frac{1}{|{x}-{y}|} k(\varepsilon x+\bar{X}_i^\varepsilon)\zeta_i^\varepsilon({x}) k(\varepsilon y+\bar{X}_i^\varepsilon)\zeta_i^\varepsilon({y})\,{\rm d}x\,{\rm d}y\\ & \quad -\mathcal{H}_l(\mathbf{x^*_1},\mathbf{x^*_2},\ldots,\mathbf{x^*_l})+o(1). \end{align*}

Since $\zeta ^\varepsilon _i\to \zeta ^*_i$ in $L^p$ weak topology and $k(\varepsilon x+\bar {X}_i^\varepsilon )\to k(\mathbf {x^*_i})$ uniformly as $\varepsilon \to 0$, we have

\begin{align*} & \iint\ln\frac{1}{|{x}-{y}|} k(\varepsilon x+\bar{X}_i^\varepsilon)\zeta_i^\varepsilon({x}) k(\varepsilon y+\bar{X}_i^\varepsilon)\zeta_i^\varepsilon({y})\,{\rm d}x\,{\rm d}y \\ & \quad =k(\mathbf{x^*_i})^2\iint\ln\frac{1}{|x-y|}\zeta^*_i(x)\zeta^*_i(y)\,{\rm d}x\,{\rm d}y+o(1), \end{align*}

where we have used the $L^p$ theory in elliptic equations and the compact embedding theorem. Thus

(3.24)\begin{align} E(\omega^\varepsilon)& =\sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^l \frac{k(\mathbf{x^*_i})^2}{4\pi}\int_{B_{R_1}(0)}\int_{B_{R_1}(0)}\ln\frac{1}{|x-y|}\zeta^*_i(x)\zeta^*_i(y)\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad -\mathcal{H}_l(\mathbf{x^*_1},\mathbf{x^*_2},\ldots,\mathbf{x^*_l})+o(1) \nonumber\\ & \leq \sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^l \frac{k(\mathbf{x^*_i})^2}{4\pi}\int_{B_{R_1}(0)}\int_{B_{R_1}(0)}\notag\\ & \quad \times \ln\frac{1}{|x-y|} \chi_{B_{\sqrt{{|d_i|}/{\pi k(\mathbf{x^*_i})}}}(0)}(x)\chi_{B_{\sqrt{{|d_i|}/{\pi k(\mathbf{x^*_i})}}}(0)}(y)\,{\rm d}x\,{\rm d}y \nonumber\\ & \quad -\mathcal{H}_l(\mathbf{x^*_1},\mathbf{x^*_2},\ldots,\mathbf{x^*_l})+o(1) \nonumber\\ & =\sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^lQ_i(k(\mathbf{x_i^*}))-\mathcal{H}_l(\mathbf{x^*_1},\mathbf{x^*_2},\ldots,\mathbf{x^*_l})+o(1) \end{align}

where the second inequality we have used (3.23) and the Riesz's rearrangement inequality.

On the other hand, for any $z_i\in B_i (i=1,\,\ldots,\,l)$, we choose a function $\hat {\omega }^\varepsilon =\sum _{i=1}^l\hat {\omega }^\varepsilon _i,\,$ where $\hat {\omega }^\varepsilon _i$ is defined by

\[ \hat{\omega}^\varepsilon_i=\frac{sgn(d_i)}{\varepsilon^2}\chi_{B_{\tau_i(\varepsilon)\varepsilon}(z_i)}. \]

Here $\tau _i(\varepsilon )\in (({1}/{2})\sqrt { {|d_i|}/{\pi k(z_i)}},\, 2\sqrt {{|d_i|}/{\pi k(z_i)}})$ is chosen to satisfy $\int _{B_i}\hat {\omega }^\varepsilon _i(x) {\rm d}\theta (x)=d_i$. Then direct calculation shows that such $\tau _i(\varepsilon )$ exists for $\varepsilon$ sufficiently small and $\lim _{\varepsilon \to 0^+}\tau _i(\varepsilon )=\sqrt {{|d_i|}/{\pi k(z_i)}}$. By the definition of $\hat {\omega }^\varepsilon$, we obtain $\hat {\omega }^\varepsilon \in \mathcal {M}_\varepsilon (D)$.

For $E(\hat {\omega }^\varepsilon )$, similarly as calculations in (3.24), we obtain

(3.25)\begin{align} E(\hat{\omega}^\varepsilon) & ={-}\frac{1}{4\pi}\sum_{i=1}^l\iint\ln|{x}-{y}| \hat{\omega}_i^\varepsilon({x}) \hat{\omega}_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y)\notag\\ & \quad - \frac{1}{2}\sum_{i=1}^l\iint h({x},{y})\hat{\omega}_i^\varepsilon({x})\hat{\omega}_i^\varepsilon({y})\theta(x)\,{\rm d}\theta(y) \nonumber\\ & \quad +\frac{1}{2}\sum_{1\leq i\neq j\leq l}\iint G(x,y)\hat{\omega}_i^\varepsilon({x})\hat{\omega}_j^\varepsilon({y})\theta(x)\,{\rm d}\theta(y) \nonumber\\ & = \sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^lQ_i(k(z_i))-\mathcal{H}_l(z_1,\ldots,z_l)+o(1). \end{align}

By $E(\omega ^\varepsilon )\geq E(\hat {\omega }^\varepsilon )$, (3.24) and (3.25), we get

\[ \sum_{i=1}^lQ_i(k(\mathbf{x^*_i}))-\mathcal{H}_l(\mathbf{x^*_1},\mathbf{x^*_2},\ldots,\mathbf{x^*_l})\geq \sum_{i=1}^lQ_i(k(z_i))-\mathcal{H}_l(z_1,\ldots,z_l)\quad \forall z_i\in B_i. \]

Thus we get (3.21). Notice that $\mathcal {H}_l(x_1,\,x_2,\,\ldots,\,x_l)-\sum _{i=1}^lQ_i(k(x_i))=\Gamma _l(x_1,\,x_2, \ldots,\,x_l)+C_0^*$, where $C_0^*$ is a universal constant. Then by the assumption that $(x_{0,1},\,\ldots,\,x_{0,l})$ is a strict local minimizer of $\Gamma _l$, we get $\mathbf {x^*_i}=x_{0,i}.$

It suffices to prove (3.22). Indeed by (3.24) we have

\begin{align*} & \iint\ln\frac{1}{|x-y|}\zeta^*_i(x)\zeta^*_i(y)\,{\rm d}x\,{\rm d}y\notag\\ & \quad = \iint\ln\frac{1}{|x-y|}\chi_{B_{\sqrt{{|d_i|}/{\pi k(\mathbf{x^*_i})}}}(0)}(x)\chi_{B_{\sqrt{{|d_i|}/{\pi k(\mathbf{x^*_i})}}}(0)}(y)\,{\rm d}x\,{\rm d}y. \end{align*}

Using strict Rearrangement inequality (see theorem 3.9, [26]), there exists a translation $\bar {\mathcal {T}}$ such that $\mathcal {T}(\zeta ^*_i)=\chi _{B_{\sqrt {{|d_i|}/{\pi k(\mathbf {x^*_i})}}}(0)}$. Notice that both the centre of $\zeta ^*_i$ and the centre of $B_{\sqrt {{|d_i|}/{\pi k(\mathbf {x^*_i})}}}(0)$ are the origin, we get $\bar {\mathcal {T}}=id$, namely, $\zeta ^*_i=\chi _{B_{\sqrt {{|d_i|}/{\pi k(\mathbf {x^*_i})}}}(0)}$.

Finally, by (3.19) we have $\int _{B_{R_1}(0)}k(\varepsilon x+\bar {X}_i^\varepsilon )(\zeta ^\varepsilon _i(x))^p\,{\rm d}x=|d_i|,\,$ which implies that

\[ \lim_{\varepsilon \to 0^+}\|\zeta^\varepsilon_i\|_{L^p}= \left(\frac{|d_i|}{k(\mathbf{x^*_i})}\right)^{{1}/{p}}=\|\zeta^*_i\|_{L^p}. \]

Using the strict convexity of $L^p$ norm, we finish the proof.

Remark 3.8 By proposition 3.7, we know that $E(\omega ^\varepsilon )$ has the following expansion

\[ E(\omega^\varepsilon) =\sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^lQ_i(k(x_{0,i}))-\mathcal{H}_l(x_{0,1},x_{0,2},\ldots,x_{0,l})+o(1). \]

Direct consequence of lemma 3.6 and proposition 3.7 is that the support set of $\omega _i^\varepsilon$ is contained in $B_i$ for $\varepsilon$ sufficiently small.

Corollary 3.9 For $\varepsilon$ sufficiently small, there holds

\[ supp(\omega_i^{\varepsilon})\subset{\subset} B_i,\quad \text{for any}\ i=1,\ldots,l. \]

Moreover, by proposition 3.7, we can repeat the classical result in [37] to show the boundary of $supp(\zeta _i^{\varepsilon })$ is a $C^1$ curve and converges to the boundary of $supp(\zeta _i^*)$ (which is a circle) in $C^1$ sense as $\varepsilon \to 0^+$, see lemma 4.10 for a detailed proof.

As a corollary of lemmas 3.4, 3.5, 3.6 and proposition 3.7, one can get the order of the functional $E_i(\omega ^\varepsilon )$ and constants $\mu _i^\varepsilon.$

Lemma 3.10 For $\varepsilon$ sufficiently small, there holds

(3.26)\begin{align} E_i(\omega^\varepsilon)& = \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+Q_i(k(x_{0,i}))-\frac{d_i^2}{2}h(x_{0,i},x_{0,i})+o(1), \end{align}

(3.27)\begin{align} \mu_i^{\varepsilon}& = \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+\frac{k(x_{0,i})}{2\pi}\int_{B_{\sqrt{{|d_i|}/{\pi k(x_{0,i})}}}(0)}\ln\frac{1}{|x^*-y|}\,{\rm d}y-|d_i|h(x_{0,i},x_{0,i})\nonumber\\ & \quad+sgn(d_i)\sum_{j=1,j\neq i}^ld_jG(x_{0,i},x_{0,j})+o(1), \end{align}

where $x^*$ is any point of $\partial B_{{\sqrt {{|d_i|}/{\pi k(x_{0,i})}}}}(0)$.

Proof. Using (3.24), one can immediately get (3.26). For (3.27), notice that for any $x\in \partial supp(\omega _i^\varepsilon )$, $\mu _i^\varepsilon =sgn(d_i)\int _{D}G(x,\,y) \omega ^\varepsilon (y) \,{\rm d}\theta (y)$. So by lemma 3.6 and proposition 3.7,

(3.28)\begin{align} \mu_i^\varepsilon& =sgn(d_i)\int_{D}G(x,y) \omega^\varepsilon(y)\,{\rm d}\theta(y) \nonumber\\ & = \frac{1}{2\pi}\int_{B_i}\ln\frac{1}{|x-y|}|\omega_i^\varepsilon(y)|\,{\rm d}\theta(y)-\int_{B_i}h(x,y)|\omega_i^\varepsilon(y)|\,{\rm d}\theta(y)\notag\\ & \quad +sgn(d_i)\sum_{j\neq i}d_jG(x_{0,i},x_{0,j})+o(1)\nonumber\\ & = \frac{1}{2\pi}\int_{B_i}\ln\frac{1}{|x-y|}|\omega_i^\varepsilon(y)|\,{\rm d}\theta(y)-|d_i|h(x_{0,i},x_{0,i})\notag\\ & \quad +sgn(d_i)\sum_{j\neq i}d_jG(x_{0,i},x_{0,j})+o(1). \end{align}

Let $x=\varepsilon x'+\bar {X}^\varepsilon _i$, then $x'\in supp(\zeta _i^\varepsilon )$. By the definition of $\zeta _i^\varepsilon$, we get

(3.29)\begin{align} & \frac{1}{2\pi}\int_{B_i}\ln\frac{1}{|x-y|}|\omega_i^\varepsilon(y)|\,{\rm d}\theta(y)\notag\\ & \quad =\frac{1}{2\pi}\int_{B_{R_1}(0)}\ln\frac{1}{\varepsilon|x'-y'|}k(\varepsilon y'+\bar{X}^\varepsilon_i)\zeta_i^\varepsilon(y')\,{\rm d}y' \nonumber\\ & \quad = \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+\frac{1}{2\pi}\int_{B_{R_1}(0)}\ln\frac{1}{|x'-y'|}k(\varepsilon y'+\bar{X}^\varepsilon_i)\zeta_i^\varepsilon(y')\,{\rm d}y'. \end{align}

By proposition 3.7 and the continuity of $k$, we have

(3.30)\begin{align} & \frac{1}{2\pi}\int_{B_{R_1}(0)}\ln\frac{1}{|x'-y'|}k(\varepsilon y'+\bar{X}^\varepsilon_i)\zeta_i^\varepsilon(y')\,{\rm d}y'\notag\\ & \quad =\frac{k(x_{0,i})}{2\pi} \int_{B_{\sqrt{{|d_i|}/{\pi k(x_{0,i})}}}(0)}\ln\frac{1}{|x^*-y'|}\,{\rm d}y'+o(1), \end{align}

where $x^*$ is any point of $\partial B_{{\sqrt {{|d_i|}/{\pi k(x_{0,i})}}}}(0)$. Taking (3.29), (3.30) into (3.28), we get (3.27).

3.3. Proof of theorem 2.3

Proof. By proposition 3.2, we know that $\omega ^\varepsilon$ has the form

\[ \omega^\varepsilon= \sum_{i=1}^l\frac{sgn(d_i)}{ \varepsilon^2} \chi_{\left\{\mathcal{G}(k\omega^\varepsilon)-\mu_i^\varepsilon>0\right\}\cap B_i}. \]

By lemma 3.6, we have $diam(supp(\omega _i^\varepsilon ))=O(\varepsilon )$. Moreover, by lemmas 3.6, 3.7 and the assumption that $(x_{0,1},\,\ldots,\,x_{0,l})$ is the strict local minimizer of $\Gamma _l$, the support set of $\omega _i^\varepsilon$ tends to $x_{0,i}$ as $\varepsilon \to 0^+$, namely,

\[ \lim_{\varepsilon \to 0^+}\sup_{x\in supp(\omega_i^\varepsilon)}|x- x_{0,i}|=0. \]

By proposition 3.7 and lemma 3.10, we get (3)(4) in theorem 2.3. The proof of theorem 2.3 is thus complete.

3.4. Proof of theorems 1.1 and 2.1

Proof. Let $l=1$ and $d_1=I$ in theorem 2.3, we get theorem 2.1. Let $\lambda ={1}/{\varepsilon ^2}$ and $u^\lambda =u^\varepsilon,\, c^\lambda =-\mu ^\varepsilon$, we get theorem 1.1.

4. Proof of theorem 2.4

Since proof of theorem 2.4 is similar to that of theorem 2.3 , we only emphasize the differences here, see proposition 4.7, lemmas 4.9, 4.10, 4.12 and 4.13.

By assumption $(f_2)$, we know that $\lim _{s\to +\infty }f(s)=+\infty$. Moreover, direct computation shows that $(f_2)$ is equivalent to

$(f_2)'.$ there exists $\delta _1\in (0,\,1)$ such that

\[ F(s) \ge \delta_1 s\, f^{{-}1}(s) \]

for any $s\geq 0$. Here $f^{-1}(s)=0$ if $t<0$ and $f^{-1}(s)$ be the inverse function of $f$ if $t\geq 0$. Let $F(s)=\int _0^s f^{-1}(t)\,{\rm d}t$.

Notice that $f^{-1}$ is nonnegative increasing continuous and $F$ is a convex $C^1$ function.

Define $F_i(s)=\int _0^sf^{-1}_i(t)\,{\rm d}t$ $(i=1,\,\ldots,\,l).$ Our idea is to consider the maximization problem

\[ (\mathcal{P}^*)\quad\mathcal{E}(\omega^\varepsilon)=\max_{\omega\in \mathcal{N}_{\varepsilon,T}(D)}\mathcal{E}(\omega), \]

where

(4.1)\begin{align} \mathcal{E}(\omega)= \frac{1}{2}\int_{D}\int_{D}G(x,y)\omega(x)\omega(y) \,{\rm d}\theta(x)\,{\rm d}\theta(y) -\frac{1}{\varepsilon^2}\sum_{i=1}^l\int_{D}F_i(sgn(d_i)\varepsilon^2\omega_i)\,{\rm d}\theta(x) \end{align}

and the set

(4.2)\begin{align} \mathcal{N}_{\varepsilon, T}(D)& =\left\{\omega=\sum_{k=1}^{l}\omega_k\in L^\infty(D)\mid \omega_k=\omega \chi_{B_k},\ \ 0\le sgn(d_i)\omega_i \le \frac{T}{ \varepsilon^2},\right. \nonumber\\ & \quad \left.\int_{D}\omega_i\,{\rm d}\theta(x) =d_i,\ \ \ i=1,\ldots,m\right\}. \end{align}

Here $T>1$ is a constant to be determined later. Note that the only difference between $\mathcal {N}_{\varepsilon, T}(D)$ and $\mathcal {M}_{\varepsilon }(D)$ defined by (3.2) is the presence of parameter $T$. However, we will show that this is a technical trick and it will not affect the final results.

4.1. Variational problem

Similarly as proof of propositions 3.1 and 3.2, we first get the existence and profile of maximizers of the functional $\mathcal {E}(\omega )$ over $\mathcal {N}_{\varepsilon, T}(D)$.

Lemma 4.1 There exists $\omega ^{\varepsilon,T}\in \mathcal {N}_{\varepsilon, T}(D)$, such that $\mathcal {E}(\omega ^{\varepsilon, T})=\sup _{\tilde {\omega }\in \mathcal {N}_{\varepsilon, T}(D)}\mathcal {E}(\tilde {\omega })$.

Proof. The proof is similar to that of proposition 3.1. So we omit it here.

Then we can get the profile of a maximizer $\omega ^{\varepsilon,T}$ as follows.

Lemma 4.2 Let $\omega ^{\varepsilon, T}$ be a maximizer defined as in lemma 4.1. Then

(4.3)\begin{align} \omega^{\varepsilon,T}=\sum_{i=1}^lsgn(d_i)\left( \frac{1}{ \varepsilon^2}f_i(\psi_i^{\varepsilon,T})\chi_{\{ 0<\psi_i^{\varepsilon,T}< f_i^{{-}1}(T) \}\cap B_i }+\frac{T}{\varepsilon^2}\chi_{\{ \psi_i^{\varepsilon,T}\geq f_i^{{-}1}(T) \}\cap B_i}\right), \end{align}

where $\psi _i^{\varepsilon,T}:=sgn(d_i)\mathcal {G}(k\omega ^{\varepsilon,T})-\sigma ^{\varepsilon,T}_i$, and $\sigma ^{\varepsilon,T}_i$ are Lagrange multipliers dependent on $\varepsilon$ for $i=1,\,\ldots,\,l$. Moreover, for $\varepsilon$ sufficiently small there holds

(4.4)\begin{equation} \sigma^{\varepsilon,T}_i\geq{-}f_i^{{-}1}(T)-C_0, \end{equation}

where $C_0>0$ is some constant independent of $\varepsilon,\, T$.

Proof. For each ${\omega }\in \mathcal {N}_{\varepsilon,T}$, we choose test functions

\[ \omega_{(s)}=\omega^{\varepsilon,T}+s({\omega}-\omega^{\varepsilon,T}),\quad \text{for} \ s\in[0,1]. \]

Since $\omega ^{\varepsilon,T}$ is a maximizer, we get $\mathcal {E}(\omega _{(s)})\leq \mathcal {E}(\omega _{(0)})$, which implies that $({{\rm d}\mathcal {E}(\omega _{(s)})})/{{\rm d}s}|_{s=0^+}\leq 0$, that is,

\begin{align*} & \int_{D}\omega \left( \mathcal{G}(k\omega^{\varepsilon,T})-\sum_{i=1}^lf_i^{{-}1}(sgn(d_i)\varepsilon^2\omega_i^{\varepsilon,T})sgn(d_i)\right)\,{\rm d}\theta(x)\\ & \quad \leq \int_{D}\omega^{\varepsilon,T}\left( \mathcal{G}(k\omega^{\varepsilon,T})-\sum_{i=1}^lf_i^{{-}1}(sgn(d_i)\varepsilon^2\omega_i^{\varepsilon,T})sgn(d_i)\right)\,{\rm d}\theta(x) \end{align*}

for all ${\omega }\in \mathcal {N}_{\varepsilon,T}.$ Using the bathtub principle, we obtain

(4.5)\begin{align} \begin{aligned} sgn(d_i) \mathcal{G}(k\omega^{\varepsilon,T}) -f_i^{{-}1}(sgn(d_i)\varepsilon^2 \omega_i^{\varepsilon,T}) & \ge \sigma_i^{\varepsilon,T} \quad \mbox{on} \left\{sgn(d_i)\omega_i^{\varepsilon,T}=\frac{T}{\varepsilon^2}\right\}\cap B_i, \\ sgn(d_i) \mathcal{G}(k\omega^{\varepsilon,T}) -f_i^{{-}1}(sgn(d_i)\varepsilon^2 \omega_i^{\varepsilon,T}) & = \sigma_i^{\varepsilon,T} \quad \mbox{on} \left\{0< sgn(d_i)\omega_i^{\varepsilon,T}<\frac{T}{\varepsilon^2}\right\}\cap B_i, \\ sgn(d_i) \mathcal{G}(k\omega^{\varepsilon,T}) -f_i^{{-}1}(sgn(d_i)\varepsilon^2 \omega_i^{\varepsilon,T}) & \le \sigma_i^{\varepsilon,T} \quad \mbox{on} \{sgn(d_i)\omega_i^{\varepsilon,T}=0\}\cap B_i, \end{aligned} \end{align}

where $\sigma _i^{\varepsilon,T}$ is a constant satisfying

(4.6)\begin{align} \sigma_i^{\varepsilon,T}=\inf\left\{s\in\mathbb{R}\mid|\{x\in B_i\mid sgn(d_i) \mathcal{G}(k\omega^{\varepsilon,T}) -f_i^{{-}1}(sgn(d_i)\varepsilon^2 \omega_i^{\varepsilon,T}) >s\}|_{\theta}\le \frac{|d_i|\varepsilon^2}{T}\right\}. \end{align}

Define $\psi _i^{\varepsilon,T}=sgn(d_i)\mathcal {G}(k\omega ^{\varepsilon,T})-\sigma ^{\varepsilon,T}_i$, then by (4.5) one has

(4.7)\begin{equation} \omega_i^{\varepsilon, T}= sgn(d_i)\left( \frac{1}{ \varepsilon^2}f_i(\psi_i^{\varepsilon,T})\chi_{\{ 0<\psi_i^{\varepsilon,T}< f_i^{{-}1}(T) \}\cap B_i }+\frac{T}{\varepsilon^2}\chi_{\{ \psi_i^{\varepsilon,T}\geq f_i^{{-}1}(T) \}\cap B_i}\right). \end{equation}

So we get (4.3).

It remains to prove that $\sigma _i^{\varepsilon, T}\geq -f_i^{-1}(T)-C_0.$ For any $x\in B_i$,

(4.8)\begin{align} & sgn(d_i) \mathcal{G}(k\omega^{\varepsilon,T}) (x) -f_i^{{-}1}(sgn(d_i)\varepsilon^2 \omega_i^{\varepsilon,T})(x)\notag\\ & \quad \geq{-}\sum_{j\neq i}|d_j|\max_{x\in B_i, y\in B_j}G(x,y) -f_i^{{-}1}(T). \end{align}

Choose $C_0=\sum _{j\neq i}|d_j|\max _{x\in B_i, y\in B_j}G(x,\,y)$. Combining (4.6) and (4.8), we get (4.4).

4.2. Asymptotic analysis of $\omega ^{\varepsilon,T}$

For simplicity, we define functionals of $\omega \in \mathcal {N}_{\varepsilon,T}$

\[ \mathcal{E}_i(\omega):= \frac{1}{2}\int_{D}G(x,y)\omega_i(x)\omega_i(y) \,{\rm d}\theta(x)\,{\rm d}\theta(y)-\frac{1}{\varepsilon^2}\int_{D}F_i(sgn(d_i)\varepsilon^2\omega_i)\,{\rm d}x. \]

Direct calculation shows that

(4.9)\begin{equation} \mathcal{E}(\omega )=\sum_{i=1}^l\mathcal{E}_i(\omega )+O(1) \end{equation}

for any $\omega \in \mathcal {N}_{\varepsilon,T}$. Here $O(1)$ is uniformly bounded about $\varepsilon$ and $T$. First we give a rough lower bound of $\mathcal {E}_i(\omega ^{\varepsilon, T})$.

Lemma 4.3 Let $\omega ^{\varepsilon,T}$ be a maximizer. Then for $i=1,\,\ldots,\,l$

(4.10)\begin{equation} \mathcal{E}_i(\omega^{\varepsilon, T})\geq \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+O(1). \end{equation}

Proof. We choose a function $\tilde {\omega }^{\varepsilon, T}=({sgn(d_i)}/{\varepsilon ^2})\chi _{B_{t(\varepsilon )\varepsilon }(x_{0,i})}+ \sum _{j\neq i}\omega _j^{\varepsilon,T}= \sum _{k=1}^l\tilde {\omega }^{\varepsilon, T}_k$. Here $t(\varepsilon )\in (({1}/{2})\sqrt { {|d_i|}/{\pi k(x_{0,i})}},\, 2\sqrt {{|d_i|}/{\pi k(x_{0,i})}})$ satisfies $\int _{B_i}\tilde {\omega }_i^{\varepsilon,T} {\rm d}\theta (x)=d_i$. Direct calculations shows that $t(\varepsilon )$ exists for $\varepsilon$ sufficiently small and $\lim _{\varepsilon \to 0^+}t(\varepsilon )=\sqrt {{|d_i|}/{\pi k(x_{0,i})}}$. Then $\tilde {\omega }^{\varepsilon, T}\in \mathcal {N}_{\varepsilon, T}(D)$. Notice that

\begin{align*} & \frac{1}{\varepsilon^2} \int_{D}F_i(sgn(d_i)\varepsilon^2\tilde{\omega}^{\varepsilon, T}_i)\,{\rm d}\theta(x)\leq \frac{1}{\varepsilon^2}F_i(1)\int_{\chi_{B_{t(\varepsilon)\varepsilon }(x_{0,i})}}\,{\rm d}\theta(x)=F_i(1)|d_i|\\ & \frac{1}{2}\int_{D}G(x,y)\tilde{\omega}^{\varepsilon, T}_i(x)\tilde{\omega}^{\varepsilon, T}_i(y) \,{\rm d}\theta(x)\,{\rm d}\theta(y)\geq \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon} +O(1), \end{align*}

so we get $\mathcal {E}(\tilde {\omega }^{\varepsilon, T})\geq ({d_i^2}/{4\pi })\ln ({1}/{\varepsilon })+\sum _{j\neq i}\mathcal {E}_j(\omega ^{\varepsilon,T})+O(1)$. By $\mathcal {E}(\omega ^{\varepsilon, T})\geq \mathcal {E}(\tilde {\omega }^{\varepsilon, T})$ and (4.9), we get (4.10).

Then we give the lower bound of Lagrange multiplier $\sigma _i^{\varepsilon, T}$.

Lemma 4.4 Let $\omega ^{\varepsilon,T}$ be a maximizer. Then for $\varepsilon$ sufficiently small there holds

(4.11)\begin{equation} \sigma_i^{\varepsilon,T}\geq{-}\frac{|d_i|}{2\pi}\ln\varepsilon-|1-2\delta_1|f_1^{{-}1}(T)-O(1). \end{equation}

Proof. By the definition of $\mathcal {E}_i$ and assumption $(f_2)'$, we get

(4.12)\begin{align} 2\mathcal{E}_i(\omega^{\varepsilon,T})& =\int_{D}\omega_i^{\varepsilon,T} \mathcal{G}(k\omega_i^{\varepsilon,T}) \,{\rm d}\theta(x) -\frac{2}{\varepsilon^2}\int_{D}F_i(sgn(d_i)\varepsilon^2\omega_i^{\varepsilon, T})\,{\rm d}\theta(x) \nonumber\\ & \leq \int_{D}\omega_i^{\varepsilon,T}( \mathcal{G}(k\omega^{\varepsilon, T})-sgn(d_i)\sigma_i^{\varepsilon,T} ) \,{\rm d}\theta(x)\nonumber\\ & \quad -2\delta_1 \int_{D}|\omega_i^{\varepsilon, T}|f_i^{{-}1}(sgn(d_i)\varepsilon^2\omega_i^{\varepsilon, T})\,{\rm d}\theta(x) +|d_i|\sigma_i^{\varepsilon, T}\nonumber\\ & = \int_{\{|\omega_i^{\varepsilon, T}|=\frac{T}{\varepsilon^2}\}}|\omega_i^{\varepsilon, T}|\left(sgn(d_i)\mathcal{G}(k\omega^{\varepsilon, T})-\sigma_i^{\varepsilon, T}-f_i^{{-}1}(T) \right) \,{\rm d}\theta(x)\nonumber\\ & \quad +(1-2\delta_1)\int_{D}|\omega_i^{\varepsilon, T}|f_i^{{-}1}( \varepsilon^2|\omega_i^{\varepsilon, T}|)\,{\rm d}\theta(x)+|d_i|\sigma_i^{\varepsilon, T} \nonumber\\ & \leq\int_{\{|\omega_i^{\varepsilon, T}|=\frac{T}{\varepsilon^2}\}}|\omega_i^{\varepsilon, T}|\left(\psi_i^{\varepsilon,T}-f_i^{{-}1}(T)-C_0 \right)_+ \,{\rm d}\theta(x)\nonumber\\ & \quad +|d_i||1-2\delta_1|f_i^{{-}1}(T)+|d_i|\sigma_i^{\varepsilon, T}+C_0|d_i|, \end{align}

where $\psi _i^{\varepsilon,T}=sgn(d_i)\mathcal {G}(k\omega ^{\varepsilon, T})-\sigma _i^{\varepsilon, T}$.

To estimate the first term of the last line of (4.12), we define $W_i^{\varepsilon,T}=(\psi _i^{\varepsilon, T}-f_i^{-1}(T) -C_0 )_+$ and $\bar {W}_i^{\varepsilon,T}=(sgn(d_i)\mathcal {G}(k\omega _i^{\varepsilon, T})-\sigma _i^{\varepsilon,T}-f_i^{-1}(T)-C_0 )_+ .$ Using (4.4), we have $\bar {W}_i^{\varepsilon,T},\, W_i^{\varepsilon,T}\in H^1_0(D)$.

On the one hand,

(4.13)\begin{equation} \int_{D} |\omega_i^{\varepsilon,T}|\bar{W}_i^{\varepsilon, T}\,{\rm d}\theta(x)=\int_{D}|\nabla\bar{W}_i^{\varepsilon,T}|^2\,{\rm d}x. \end{equation}

On the other hand, by (4.3), we get $supp(W_i^{\varepsilon, T})\cap B_i\subseteq \{|\omega _i^{\varepsilon,T}|={T}/{\varepsilon ^2}\}$, which implies that

\begin{align*} \int_{D}|\omega_i^{\varepsilon,T}|\bar{W}_i^{\varepsilon,T}\,{\rm d}\theta(x) & \leq \int_{D}|\omega_i^{\varepsilon,T}|W_i^{\varepsilon,T}\,{\rm d}\theta(x)+O(1) \nonumber\\ & =\int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}|\omega_i^{\varepsilon,T}|W_i^{\varepsilon, T}\,{\rm d}\theta(x)+O(1). \end{align*}

Using the Sobolev imbedding $W^{1,1}(B_i)\subset L^2(B_i)$, we get

\begin{align*} & \int_{D}|\omega_i^{\varepsilon, T}|\bar{W}_i^{\varepsilon,T}\,{\rm d}\theta(x) \nonumber\\ & \quad \leq\frac{T}{\varepsilon^2}\int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}W_i^{\varepsilon,T}\,{\rm d}\theta(x)+O(1)\\ & \quad \leq \frac{CT}{\varepsilon^2}\left|\left\{|\omega_i^{\varepsilon, T}|=\frac{T}{\varepsilon^2}\right\}\right|^{{1}/{2}}\left( \int_{B_i}W_i^{\varepsilon, T}\,{\rm d}x+\int_{B_i}|\nabla W_i^{\varepsilon,T}|\,{\rm d}x\right)+ O(1)\\ & \quad \leq C \left|\left\{|\omega_i^{\varepsilon, T}|=\frac{T}{\varepsilon^2}\right\}\right|^{{1}/{2}} \int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}sgn(d_i)\omega_i^{\varepsilon, T}\bar{W}_i^{\varepsilon, T}\,{\rm d}x \\ & \quad \quad +\frac{CT}{\varepsilon^2}\left|\left\{|\omega_i^{\varepsilon, T}|=\frac{T}{\varepsilon^2}\right\}\right|^{{1}/{2}}\int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}|\nabla\bar{W}_i^{\varepsilon,T}|\,{\rm d}x+O(1). \end{align*}

So by the fact that $|\{|\omega _i^{\varepsilon,T}|={\Lambda }/{\varepsilon ^2}\}|=O(\varepsilon ^2)$ and Hölder's inequality, for $\varepsilon$ sufficiently small we get

(4.14)\begin{align} & \int_{D}|\omega_i^{\varepsilon, T}|\bar{W}_i^{\varepsilon,T}\,{\rm d}\theta(x) \nonumber\\ & \quad \leq\frac{CT}{\varepsilon^2}\left|\left\{|\omega_i^{\varepsilon, T}|=\frac{T}{\varepsilon^2}\right\}\right|^{{1}/{2}}\int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}|\nabla\bar{W}_i^{\varepsilon,T}|\,{\rm d}x+O(1) \nonumber\\ & \quad \leq\frac{CT}{\varepsilon^2}\left|\left\{|\omega_i^{\varepsilon, T}|=\frac{T}{\varepsilon^2}\right\}\right|\left( \int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}|\nabla\bar{W}_i^{\varepsilon, T}|^2\,{\rm d}x\right)^{{1}/{2}} +O(1) \nonumber\\ & \quad \leq C|d_i|\left( \int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}|\nabla\bar{W}_i^{\varepsilon, T}|^2\,{\rm d}x\right)^{{1}/{2}} +O(1). \end{align}

Combining (4.13) and (4.14), we get $\int _{D}|\omega _i^{\varepsilon, T}|\bar {W}_i^{\varepsilon,T}\,{\rm d}\theta (x)=O(1),\,$ which implies that

(4.15)\begin{equation} \int_{\{|\omega_i^{\varepsilon, T}|={T}/{\varepsilon^2}\}}|\omega_i^{\varepsilon, T}|W_i^{\varepsilon,T}\,{\rm d}\theta(x)=O(1). \end{equation}

Taking (4.15) into (4.12) and using lemma 4.3, we get (4.11).

Using assumptions $(f_1)$–$(f_3)$ and the rearrangement inequality, one can get the upper bound of $\psi _i^{\varepsilon, T}$. As a result, we show that the vortex patch part of $\omega _i^{\varepsilon,T}$ indeed vanishes by choosing $T$ sufficiently large.

Lemma 4.5 Let $\omega ^{\varepsilon,T}$ be a maximizer as in lemma 4.1. Then

(4.16)\begin{equation} \|\psi_i^{\varepsilon,T}\|_{L^\infty(B_i)}\leq |1-2\delta_1|f_i^{{-}1}(T)+\frac{|d_i|}{4\pi}\ln T+O(1). \end{equation}

As a consequence, one can choose $T=T_0$ sufficiently large such that

\[ |\{ \psi_i^{\varepsilon, T_0}\geq f_i^{{-}1}(T_0) \}\cap B_i|=0, \]

and so $\omega _i^{\varepsilon, T_0}$ has the form

(4.17)\begin{equation} \omega_i^{\varepsilon, T_0}=sgn(d_i)\frac{1}{ \varepsilon^2}f_i(\psi_i^{\varepsilon,T_0})\chi_{\{ \psi_i^{\varepsilon,T_0}>0 \}\cap B_i }. \end{equation}

Proof. For any $x\in B_i$, using the definition of $\psi _i^{\varepsilon, T}$ and the rearrangement inequality

\begin{align*} \psi_i^{\varepsilon, T}(x)& =sgn(d_i)\mathcal{G}(k\omega^{\varepsilon,T})(x)-\sigma^{\varepsilon,T}_i \\ & \leq\frac{1}{2\pi}\int_{D}\ln\frac{1}{|x-y|}|\omega_i^{\varepsilon, T}|(y)\,{\rm d}y\theta(y)-\sigma_i^{\varepsilon, T}+O(1)\\ & \leq \frac{k_1T}{2\pi\varepsilon^2}\int_{B_{\sqrt{{|d_i|}/{(\pi k_1T)} }\varepsilon}(0)}\ln\frac{1}{|y|}\,{\rm d}y-\sigma_i^{\varepsilon,T}+O(1)\\ & = \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+\frac{|d_i|}{4\pi}\ln T-\sigma_i^{\varepsilon,T}+O(1). \end{align*}

Thus using lemma 4.4, we obtain

\[ \psi_i^{\varepsilon, \Lambda}(x)\leq |1-2\delta_1|f_i^{{-}1}(T)+\frac{|d_i|}{4\pi}\ln T+O(1). \]

It follows from assumption $(f_3)$ that for each $a_0 > 0$, $\lim _{s\to +\infty }f_i(s)e^{-a_0 s}=0,\,$ which implies that $\lim _{s\to +\infty } \tau f^{-1}(s)-\ln s=+\infty.$ Thus we can choose $T=T_0$ sufficiently large such that

\[ (1-|1-2\delta_1|)f_i^{{-}1}(T_0)>\frac{|d_i|}{4\pi}\ln T_0+O(1), \]

that is, $|1-2\delta _1|f_i^{-1}(T_0)+({|d_i|}/{4\pi })\ln T_0+O(1)< f_i^{-1}(T_0)$. Thus we have $|\{\psi _i^{\varepsilon,T}\geq f_i^{-1}(T) \}\cap B_i|=0.$ Using lemma 4.2, we get (4.17).

In the following, we shall abbreviate $(\mathcal {N}_{\varepsilon, T_0}(D);\omega _i^{\varepsilon, T_0};\sigma _i^{\varepsilon, T_0};\psi _i^{\varepsilon,T_0})$ as $(\mathcal {N}_{\varepsilon }(D);\omega _i^{\varepsilon };\sigma _i^{\varepsilon };\psi _i^{\varepsilon })$ for $i=1,\,\ldots,\,l$. By lemma 4.5, we know that any maximizer $\omega ^{\varepsilon }$ of the maximization problem $(\mathcal {P}^*)$ has the form of (2.5).

Similarly as lemma 3.6, we can get the diameter of $supp(\omega _i^\varepsilon )$ is of the order $O(\varepsilon ).$

Lemma 4.6 Let $\omega ^{\varepsilon }$ be a maximizer. Then for $\varepsilon$ sufficiently small, there holds

(4.18)\begin{equation} \bar{r}_1\varepsilon\leq diam\left(supp( \omega_i^\varepsilon)\right)\le \bar{R}_1\varepsilon \end{equation}

for some $0<\bar {r}_1< \bar {R}_1$ independent of $\varepsilon$.

Proof. Since $|\omega _i^\varepsilon |\leq {T_0}/{\varepsilon ^2}$ and $\int _{B_i}\omega _i^\varepsilon \,{\rm d}\theta (x)=d_i$, we get $|supp( \omega _i^\varepsilon )|\geq C\varepsilon ^2$, which implies that $diam(supp( \omega _i^\varepsilon ))\geq \bar {r}_1\varepsilon$ for some $\bar {r}_1>0$.

Similarly as the proof of lemma 3.6, one can get the existence of $\bar {R}_1>1$ such that

\[ diam\left(supp( \omega_i^\varepsilon)\right)\leq \bar{R}_1\varepsilon.\]

□

Finally, we analyse the limiting location of $\omega _i^\varepsilon$ as $\varepsilon \to 0^+$, which is the most important part in our construction. To this end we define the centre of $\omega _i^\varepsilon$ by

\[ \hat{X}_i^\varepsilon:=\frac{1}{d_i}\int_{D}x \omega_{i}^\varepsilon(x) \,{\rm d}\theta(x)\quad \forall\ i=1,\ldots,l. \]

Since $\overline {B_i}$ is compact, we may choose a subsequence $\{\hat {X}_i^{\varepsilon _n}\}_{n=1}^{ \infty }$(still denoted by $\bar {X}_i^\varepsilon$) satisfying

\[ \lim_{\varepsilon\to 0^+}\hat{X}_i^\varepsilon=x^*_i\in \overline{B_i}. \]

Define the scaled functions

(4.19)\begin{equation} \xi^\varepsilon_i(x)=sgn(d_i)\varepsilon^2\omega_i^\varepsilon(\varepsilon x+\hat{X}_i^\varepsilon)\quad x\in D_\varepsilon. \end{equation}

Here $D_\varepsilon =\{x\in \mathbb {R}^2 \mid \varepsilon x+\hat {X}_i^\varepsilon \in D\}$.

By lemma 4.6, we know that the support set of $\xi ^\varepsilon _i$ is contained in $B_{\bar { R}_1}(0)$. Notice that

(4.20)\begin{equation} |d_i|=\int_{B_i}|\omega_i|(x)\,{\rm d}\theta(x)=\int_{B_{\bar{ R}_1}(0)}k(\varepsilon x+\hat{X}_i^\varepsilon)\xi^\varepsilon_i(x)\,{\rm d}x, \end{equation}

which implies that

(4.21)\begin{equation} \lim_{\varepsilon \to 0^+}\int_{B_{\bar{ R}_1}(0)}\xi_i^\varepsilon(x)\,{\rm d}x=\frac{|d_i|}{k(x^*_i)}. \end{equation}

Since $\|\xi ^\varepsilon _i\|_{L^\infty (B_{\bar { R}_1}(0))}\leq T_0$, $\xi ^\varepsilon _i$ is uniformly bounded in $L^p(B_{\bar { R}_1}(0))$ for any $p\in [1,\,+\infty ]$. So up to a subsequence, we may assume that $\xi ^\varepsilon _i\to \xi ^*_i$ in $L^p$ weak topology as $\varepsilon \to 0.$ By the definition of $\hat {X}_i^\varepsilon$ and (4.19), one can get

\[ \int_{B_{\bar{ R}_1}(0)}x\xi_i^*(x)\,{\rm d}x=\lim_{\varepsilon \to 0^+}\int_{B_{\bar{ R}_1}(0)}x\xi_i^\varepsilon(x)\,{\rm d}x=0. \]

Since $\omega ^\varepsilon$ is a maximizer of $\mathcal {E}$, using lemma 4.6 we can get the necessary condition of $x^*_i$.

Proposition 4.7 There holds

(4.22)\begin{equation} \mathcal{H}_l(x^*_1,\ldots,x^*_l)-\sum_{i=1}^l\frac{d_i^2}{8\pi}\ln k(x^*_i)=\min_{x_i\in B_i}\mathcal{H}_l(x_1,x_2,\ldots,x_l)-\sum_{i=1}^l\frac{d_i^2}{8\pi}\ln k(x_i). \end{equation}

As a consequence, $x^*_i =x_{0,i}$ for $i=1,\,\ldots,\,l.$

Proof. On the one hand, using lemma 4.6 and the definition of $x^*_i$, we get

\begin{align*} & \frac{1}{2}\int_{D}\int_{D}h({x},{y})\omega_i^\varepsilon({x})\omega_i^\varepsilon({y})\theta(x)\,{\rm d}\theta(y)=\frac{d_i^2}{2}h(x_i^*,x_i^*)+o(1),\\ & \sum_{1\leq i\neq j\leq l}\iint G(x,y)\omega_i^\varepsilon({x})\omega_j^\varepsilon({y})\theta(x)\,{\rm d}\theta(y)=\sum_{1\leq i\neq j\leq l}d_id_jG(x_i^*,x_j^*)+o(1). \end{align*}

By the definition of $\xi _i^\varepsilon$, we have

\begin{align*} & -\frac{1}{4\pi}\iint\ln|{x}-{y}| \omega_i^\varepsilon({x}) \omega_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y)\\ & \quad = \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\frac{1}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}\xi_i^\varepsilon({x})\xi_i^\varepsilon({y})k(\varepsilon x+\hat{X}_i^\varepsilon)k(\varepsilon y+\hat{X}_i^\varepsilon)\,{\rm d}x\,{\rm d}y. \end{align*}

Since $supp(\xi _i^\varepsilon )\subseteq B_{\bar { R}_1}(0)$ and $||\xi ^\varepsilon _i||_{L^\infty }\leq T_0$, we obtain

(4.23)\begin{align} & -\frac{1}{4\pi}\iint\ln|{x}-{y}| \omega_i^\varepsilon({x}) \omega_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y)\nonumber\\ & \quad =\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\frac{k(x_i^*)^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}\xi^\varepsilon_i(x)\xi^\varepsilon_i(y)\,{\rm d}x\,{\rm d}y+o(1). \end{align}

For the term involving $F_i$, we have

\begin{align*} \frac{1}{\varepsilon^2}\int_{D}F_i(sgn(d_i)\varepsilon^2\omega_i^\varepsilon)\,{\rm d}\theta(x)& =\int_{B_{\bar{ R}_1}(0)}F_i(\xi_i^\varepsilon)(x)k(\varepsilon x+\hat{X}_i^\varepsilon)\,{\rm d}x \nonumber\\ & =k(x_i^*)\int_{\mathbb{R}^2}F_i(\xi_i^\varepsilon)(x)\,{\rm d}x+o(1). \end{align*}

Taking those into the definition of $\mathcal {E}$, we get

(4.24)\begin{align} \mathcal{E}(\omega^\varepsilon) & = \sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^l\frac{k(x_i^*)^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}\xi^\varepsilon_i(x)\xi^\varepsilon_i(y)\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad -\mathcal{H}_l(x_1^*,\ldots,x_l^*) -\sum_{i=1}^lk(x_i^*)\int_{\mathbb{R}^2}F_i(\xi_i^\varepsilon)(x)\,{\rm d}x+o(1). \end{align}

On the other hand, for any $z_i\in B_i (i=1,\,\ldots,\,l)$, we choose test functions $\hat {\omega }^\varepsilon =\sum _{i=1}^l\hat {\omega }^\varepsilon _i,\,$ where $\hat {\omega }^\varepsilon _i$ is defined by

\[ \hat{\omega}^\varepsilon_i=\frac{sgn(d_i)}{\varepsilon^2}\xi_i^\varepsilon\left(\frac{\cdot{-}z_i}{\bar{\tau}_i(\varepsilon)\varepsilon}\right). \]

$\bar {\tau }_i(\varepsilon )\in (({1}/{2})\sqrt {{k(x_i^*)}/{k(z_i)}} ,\, 2\sqrt {{k(x_i^*)}/{k(z_i)}} )$ is chosen to satisfy $\int _{B_i}\hat {\omega }^\varepsilon _i(x)\,{\rm d}\theta (x)=d_i$. By (4.20) and (4.21), one can prove that such $\bar {\tau }_i(\varepsilon )$ exists for $\varepsilon$ sufficiently small and $\bar {\tau }_i(\varepsilon )=\sqrt {({k(x_i^*)}/{k(z_i)})}+o(1)$. Now we calculate the energy expansion of $\hat {\omega }^\varepsilon$. It is not hard to prove that

\begin{align*} & - \frac{1}{2}\sum_{i=1}^l\iint h({x},{y})\hat{\omega}_i^\varepsilon({x})\hat{\omega}_i^\varepsilon({y})\theta(x)\,{\rm d}\theta(y) \\ & \quad +\frac{1}{2}\sum_{1\leq i\neq j\leq l}\iint G(x,y)\hat{\omega}_i^\varepsilon({x})\hat{\omega}_j^\varepsilon({y})\theta(x)\,{\rm d}\theta(y)\\ & \qquad ={-}\mathcal{H}_l(z_1,\ldots,z_l)+o(1). \end{align*}

Similar as (4.23), we get

(4.25)\begin{align} & -\frac{1}{4\pi}\iint\ln|{x}-{y}| \hat{\omega}_i^\varepsilon({x}) \hat{\omega}_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y) \nonumber\\ & \quad = \frac{d_i^2}{4\pi}\ln\frac{1}{\bar{\tau}_i(\varepsilon)\varepsilon}+\frac{k(z_i)^2\bar{\tau}_i(\varepsilon)^4}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}\xi^\varepsilon_i(x)\xi^\varepsilon_i(y)\,{\rm d}x\,{\rm d}y+o(1). \end{align}

Taking $\bar {\tau }_i(\varepsilon )=\sqrt {({k(x_i^*)}/{k(z_i)})}+o(1)$ into (4.25) we obtain

\begin{align*} & -\frac{1}{4\pi}\iint\ln|{x}-{y}| \hat{\omega}_i^\varepsilon({x}) \hat{\omega}_i^\varepsilon({y})\,{\rm d}\theta(x)\,{\rm d}\theta(y)\\ & \quad = \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\frac{d_i^2}{8\pi}\ln k(z_i)-\frac{d_i^2}{8\pi}\ln k(x_i^*)\\ & \qquad +\frac{k(x_i^*)^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}\xi^\varepsilon_i(x)\xi^\varepsilon_i(y)\,{\rm d}x\,{\rm d}y+o(1). \end{align*}

For the term involving $F_i$, we get

\begin{align*} \frac{1}{\varepsilon^2}\int_{D}F_i(sgn(d_i)\varepsilon^2\hat{\omega}_i^\varepsilon)\,{\rm d}\theta(x)& =\int_{B_{\bar{ R}_1}(0)}F_i(\xi_i^\varepsilon)(x')k(\bar{\tau}_i(\varepsilon)\varepsilon x'+z_i)\bar{\tau}_i(\varepsilon)^2\,{\rm d}x'\\ & = k(x_i^*)\int_{\mathbb{R}^2}F_i(\xi_i^\varepsilon)(x)\,{\rm d}x+o(1). \end{align*}

Taking those into the definition of $\mathcal {E}(\hat {\omega }^\varepsilon )$, we get

(4.26)\begin{align} \mathcal{E}(\hat{\omega}^\varepsilon) & =\sum_{i=1}^l\left( \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\frac{d_i^2}{8\pi}\ln k(z_i)-\frac{d_i^2}{8\pi}\ln k(x_i^*)\right.\nonumber\\ & \quad \left.+\frac{k(x_i^*)^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}\xi^\varepsilon_i(x)\xi^\varepsilon_i(y)\,{\rm d}x\,{\rm d}y\right) \nonumber\\ & \quad -\mathcal{H}_l(z_1,\ldots,z_l)-\sum_{i=1}^lk(x_i^*)\int_{\mathbb{R}^2}F_i(\xi_i^\varepsilon)(x)\,{\rm d}x+o(1). \end{align}

Since $\mathcal {E}(\omega ^\varepsilon )\geq \mathcal {E}(\hat {\omega }^\varepsilon )$, by (4.24) and (4.26) we get

\[ \sum_{i=1}^l\frac{d_i^2}{8\pi}\ln k(x^*_i)-\mathcal{H}_l(x^*_1,\ldots,x^*_l)\geq \sum_{i=1}^l\frac{d_i^2}{8\pi}\ln k(z_i)-\mathcal{H}_l(z_1,\ldots,z_l)\quad \forall\ z_i\in B_i. \]

Thus we get (4.22). By the assumption that $(x_{0,1},\,\ldots,\,x_{0,l})$ is a strict local minimizer of $\Gamma _l$, we have $x^*_i=x_{0,i}.$

Remark 4.8 By lemma 4.6 and proposition 4.7, we get for $\varepsilon$ sufficiently small,

\[ supp(\omega_i^{\varepsilon})\subset{\subset} B_i, \quad \forall i=1,\ldots,l. \]

Using lemma 4.6 and proposition 4.7, we can further get the accurate estimates of $\mathcal {E}$ and $\sigma _i^\varepsilon$. To this end, we define the scaled functions of $\psi _i$ by

\[ \Psi_i^\varepsilon(x)=(sgn(d_i)\mathcal{G}(k\omega^\varepsilon)-\sigma_i^\varepsilon)(\varepsilon x+\hat{X}_i^\varepsilon)\quad x\in (B_i)_\varepsilon. \]

Here $( B_i)_\varepsilon =\{x\in \mathbb {R}^2 \mid \varepsilon x+\hat {X}_i^\varepsilon \in B_i\}$. By lemma 4.2 and remark 4.8, we have $supp(\xi ^\varepsilon _i)=supp((\Psi _i^\varepsilon )_+)\subseteq B_{\bar { R}_1}(0)$.

It follows from (4.17) that $\Psi _i^\varepsilon$ satisfies

(4.27)\begin{equation} \xi^\varepsilon_i(x)=f_i(\Psi_i^\varepsilon)(x)\quad x\in (B_i)_\varepsilon. \end{equation}

Thus by the definition of $\Psi _i^\varepsilon$ and $\xi ^\varepsilon _i$, we get

(4.28)\begin{align} -\Delta\Psi_i^\varepsilon(x)& =sgn(d_i)\varepsilon^2k(\varepsilon x+\hat{X}_i^\varepsilon)\omega_i^\varepsilon(\varepsilon x+\hat{X}_i^\varepsilon) \nonumber\\ & =k(\varepsilon x+\hat{X}_i^\varepsilon)\xi^\varepsilon_i(x)=k(\varepsilon x+\hat{X}_i^\varepsilon)f_i(\Psi_i^\varepsilon)(x). \end{align}

By (4.21) and (4.27), we have

(4.29)\begin{equation} \lim_{\varepsilon \to 0^+}\int_{B_{\bar{ R}_1}(0)}f_i(\Psi_i^\varepsilon)(x)\,{\rm d}x=\lim_{\varepsilon \to 0^+}\int_{B_{\bar{ R}_1}(0)}\xi_i^\varepsilon(x)\,{\rm d}x=\frac{|d_i|}{k(x^*_i)}. \end{equation}

Let $\Psi ^*_i$ be the unique radial function satisfying

(4.30)\begin{equation} \begin{cases} \displaystyle-\Delta\Psi_i^*(x)=k(x^*_i)f_i(\Psi_i^*)(x)\quad x\in \mathbb{R}^2,\\ \displaystyle\int_{\mathbb{R}^2}f_i(\Psi_i^*)(x)\,{\rm d}x=\frac{|d_i|}{k(x^*_i)}. \end{cases} \end{equation}

Then $\Psi ^*_i(x)=U_i(k(x^*_i)^{{1}/{2}}x)$, where $U_i$ is the unique radial function satisfying

(4.31)\begin{equation} \begin{cases} -\Delta U_i(x)= f_i(U_i)(x)\quad x\in \mathbb{R}^2,\\ \displaystyle\int_{\mathbb{R}^2}f_i(U_i)(x)\,{\rm d}x=|d_i|. \end{cases} \end{equation}

We first show that $\xi _i^*$ is a radial function. Denote $\tilde {\xi }_i^\varepsilon$ the radially symmetric decreasing Lebesque rearrangement function of $\xi _i^\varepsilon$. Up to a subsequence we may assume that $\tilde {\xi }_i^\varepsilon \to \tilde {\xi }_i^*$ weakly in $L^p(B_{\bar { R}_1}(0))$ as $\varepsilon \to 0^+$.

Lemma 4.9 There holds

\[ \xi_i^*=\tilde{\xi}^*_i. \]

So $\xi _i^*$ is a radially symmetric function.

Proof. On the one hand, by the Riesz's rearrangement inequality, we have

\begin{align*} & \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\xi_i^\varepsilon(x)\xi_i^\varepsilon(y)\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad \le \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\tilde{\xi}^\varepsilon(x)\tilde{\xi}^\varepsilon(y)\,{\rm d}x\,{\rm d}y, \end{align*}

which implies that

(4.32)\begin{align} & \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\xi_i^*(x)\xi_i^*(y)\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad \le \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\tilde{\xi}_i^*(x)\tilde{\xi}_i^*(y)\,{\rm d}x\,{\rm d}y. \end{align}

On the other hand, let $\tilde {\omega }^\varepsilon =\tilde {\omega }^\varepsilon _i+\sum _{j=1,j\neq i}^l\omega ^\varepsilon _j\in \mathcal {N}_\varepsilon (D)$ satisfying

\[ \tilde{\omega}^\varepsilon_i(x)=\left\{ \begin{array}{ll} \displaystyle\dfrac{sgn(d_i)}{\varepsilon^2}\tilde{\xi}_i^\varepsilon \left(\dfrac{x-\hat{X}_i^\varepsilon}{c(\varepsilon)\varepsilon}\right) & \text{if} \ x\in B_{\bar{ R}_1}(\hat{X}_i^\varepsilon), \\ 0 & \text{if} \ x\in D\backslash B_{\bar{ R}_1}(\hat{X}_i^\varepsilon), \end{array} \right. \]

where $c(\varepsilon )$ is a constant such that $\int _{B_i}\tilde {\omega }^\varepsilon _i\,{\rm d}\theta (x)=d_i$. Then $c(\varepsilon )=1+o(1)$. Similarly as the proof of proposition 4.7, we get

\begin{align*} \mathcal{E}({\omega}^\varepsilon) & = \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+ \frac{d_i^2}{4\pi}\int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\xi_i^\varepsilon(x)\xi_i^\varepsilon(y)\,{\rm d}x\,{\rm d}y\\ & \quad -\frac{d_i^2}{2}h(x^*_i,x^*_i)+\sum_{j=1,j\neq i}^ld_id_jG(x^*_i,x^*_j)\\ & \quad -\frac{1}{\varepsilon^2}\int_{B_i}F_i(sgn(d_i)\varepsilon^2\omega^\varepsilon_i)\,{\rm d}\theta(x)+\mathcal{E}\left(\sum_{j=1,j\neq i}^l\omega_j^\varepsilon\right)+o(1),\\ \mathcal{E}(\tilde{\omega}^\varepsilon) & = \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+ \frac{d_i^2}{4\pi}\int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\tilde{\xi}_i^\varepsilon(x)\tilde{\xi}_i^\varepsilon(y)\,{\rm d}x\,{\rm d}y\\ & \quad -\frac{d_i^2}{2}h(x^*_i,x^*_i)+\sum_{j=1,j\neq i}^ld_id_jG(x^*_i,x^*_j)\\ & \quad -\frac{1}{\varepsilon^2}\int_{B_i}F_i(sgn(d_i)\varepsilon^2\tilde{\omega}^\varepsilon_i)\,{\rm d}\theta(x)+\mathcal{E} \left(\sum_{j=1,j\neq i}^l\omega_j^\varepsilon\right)+o(1), \end{align*}

and

\[ \frac{1}{\varepsilon^2}\int_{B_i}F_i(sgn(d_i)\varepsilon^2\omega^\varepsilon_i)\,{\rm d}\theta(x)=\frac{1}{\varepsilon^2}\int_{B_i}F_i(sgn(d_i)\varepsilon^2\tilde{\omega}^\varepsilon_i)\,{\rm d}\theta(x)+o(1). \]

Since $\mathcal {E}(\tilde {\omega }^\varepsilon )\le \mathcal {E}({\omega }^\varepsilon )$, we conclude that

\begin{align*} & \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\xi_i^\varepsilon(x)\xi_i^\varepsilon(y)\,{\rm d}x\,{\rm d}y\\ & \quad \ge \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\tilde{\xi}^\varepsilon(x)\tilde{\xi}^\varepsilon(y)\,{\rm d}x\,{\rm d}y+o(1), \end{align*}

which implies that

\begin{align*} & \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\xi_i^*(x)\xi_i^*(y)\,{\rm d}x\,{\rm d}y\\ & \quad \geq \int_{B_{\bar{ R}_1}(0)}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x-y|}\tilde{\xi}_i^*(x)\tilde{\xi}_i^*(y)\,{\rm d}x\,{\rm d}y. \end{align*}

Thus the equality holds in (4.32). By the strict Riesz's rearrangement inequality (see theorem 3.9, [26]), there exists a translation $\mathcal {T}$ such that $\mathcal {T}\xi _i^*=\tilde {\xi }_i^*$. Since

\[ \int_{B_{\bar{ R}_1}(0)}x\xi_i^*(x)\,{\rm d}x=\int_{B_{\bar{ R}_1}(0)}x\tilde{\xi}_i^*(x)\,{\rm d}x=0, \]

we get $\xi _i^*=\tilde {\xi }^*_i$.

Lemma 4.10 There holds as $\varepsilon \to 0$,

\[ \Psi_i^\varepsilon\to \Psi^*_i\quad \text{in}\ C^1_{loc}(\mathbb{R}^2). \]

As a consequence, for $\varepsilon$ sufficiently small, $\partial (supp(\xi _i^\varepsilon ))=\{x\in B_{\bar { R}_1}(0)\mid \Psi _i^\varepsilon (x)=0 \}$ is a $C^1$ curve and converges to the circle $\{x\in B_{\bar { R}_1}(0)\mid \Psi _i^*(x)=0 \}$ as $\varepsilon \to 0$.

Proof. For any $R>\bar { R}_1$, notice that $\xi _i^\varepsilon$ is uniformly bounded in $L^\infty (B_{2R}(0))$. Thus, by (4.28) and classical elliptic estimates, $\Psi ^\varepsilon _i$ is uniformly bounded in $W^{2,p}(B_{R}(0))$ for every $1\le p<+\infty$. By the Sobolev embedding theorem, we may conclude that $\Psi ^\varepsilon _i$ is compact in $C^{1,\alpha }(B_{R}(0))$ for every $0<\alpha <1$. Then up to a subsequence we may assume $\Psi ^\varepsilon \to \Psi$ in $C^{1,\alpha }(B_{R}(0))$. By (4.28) and (4.29), $\Psi$ satisfies

\[ \begin{cases} -\Delta \Psi =k(x^*_i)\xi_i^*=k(x^*_i)f_i(\Psi)\quad x\in\mathbb{R}^2,\\ \displaystyle\int_{\mathbb{R}^2}f_i(\Psi) \,{\rm d}x=\frac{|d_i|}{k(x^*_i)}. \end{cases} \]

By lemma 4.9, $\xi _i^*$ is a radial function. Using the Green's function representation, $\Psi$ is also radial. By the uniqueness of radial solutions of (4.30), we have $\Psi =\Psi ^*_i$.

By the strong maximum principle, one can show that for any $x\in \{x\in B_{\bar { R}_1}(0)\mid \Psi _i^*(x)=0 \}$, $|\nabla \Psi _i^*(x)|\neq 0$. Thus by the implicit function theorem, we get that for $\varepsilon$ sufficiently small, $\{x\in B_{\bar { R}_1}(0)\mid \Psi _i^\varepsilon (x)=0 \}$ is a $C^1$ curve and converges to $\{x\in B_{\bar { R}_1}(0)\mid \Psi _i^*(x)=0 \}$ in $C^1$ sense as $\varepsilon \to 0$.

As a corollary, we have $\xi _i^\varepsilon \to \xi _i^*$ in $L^p$ topology as $\varepsilon \to 0^+$.

Corollary 4.11 For any $p>1$, there holds as $\varepsilon \to 0$

\[ \xi_i^\varepsilon\to \xi^*_i=f_i(\Psi_i^*)\quad \text{in}\ L^p_{loc}(\mathbb{R}^2). \]

Proof. Using (4.27) and lemma 4.10, one can immediately get the result.

Using lemma 4.10 and corollary 4.11, we can get the asymptotic expansion of the energy $\mathcal {E}(\omega _i^\varepsilon )$ as follows

Lemma 4.12 There holds

(4.33)\begin{equation} \mathcal{E}(\omega^\varepsilon) = \sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^l\frac{d_i^2}{8\pi}\ln k(x_{0,i})-\mathcal{H}_l(x_{0,1},\ldots,x_{0,l})+C_0+o(1), \end{equation}

where $C_0=({1}/{4\pi })\sum _{i=1}^l\int _{\mathbb {R}^2}\int _{\mathbb {R}^2}\ln ({1}/{|x-y|})f_i\circ U_i(x)f_i\circ U_i(y)\,{\rm d}x\,{\rm d}y-\sum _{i=1}^l\int _{\mathbb {R}^2} F_i\circ f_i\circ U_i(x)\,{\rm d}x$ is a constant independent of $\varepsilon$.

Proof. The proof is based on the proof of proposition 4.7. By (4.24), we have

(4.34)\begin{align} \mathcal{E}(\omega^\varepsilon) & = \sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^l\frac{k(x_{0,i})^2}{4\pi}\iint\ln\frac{1}{|x-y|}\xi^\varepsilon_i(x)\xi^\varepsilon_i(y)\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad -\mathcal{H}_l(x_{0,1},\ldots,x_{0,l}) -\sum_{i=1}^lk(x_{0,i})\int_{\mathbb{R}^2}F_i(\xi_i^\varepsilon)(x)\,{\rm d}x+o(1). \end{align}

By lemma 4.10 and corollary 4.11, we obtain

\begin{align*} & \frac{k(x_{0,i})^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}\xi^\varepsilon_i(x)\xi^\varepsilon_i(y)\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad =\frac{k(x_{0,i})^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}f_i\circ\Psi_i^*(x)f_i\circ\Psi_i^*(y)\,{\rm d}x\,{\rm d}y+o(1),\\ & k(x_{0,i})\int_{\mathbb{R}^2}F_i(\xi_i^\varepsilon)(x)\,{\rm d}x=k(x_{0,i})\int_{\mathbb{R}^2}F_i\circ f_i\circ\Psi_i^*(x)\,{\rm d}x+o(1). \end{align*}

Taking those into (4.34), we have

(4.35)\begin{align} \mathcal{E}(\omega^\varepsilon) & = \sum_{i=1}^l\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\sum_{i=1}^l\frac{k(x_{0,i})^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}f_i\circ\Psi_i^*(x)f_i\circ\Psi_i^*(y)\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad-\mathcal{H}_l(x_{0,1},\ldots,x_{0,l}) -\sum_{i=1}^lk(x_{0,i})\int_{\mathbb{R}^2}F_i\circ f_i\circ\Psi_i^*(x)\,{\rm d}x+o(1). \end{align}

However, by the definition of $\Psi _i^*$ and (4.31), we get

(4.36)\begin{align} & \sum_{i=1}^l\frac{k(x_{0,i})^2}{4\pi}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{1}{|x-y|}f_i\circ\Psi_i^*(x)f_i\circ\Psi_i^*(y)\,{\rm d}x\,{\rm d}y\nonumber\\ & \quad-\sum_{i=1}^lk(x_{0,i})\int_{\mathbb{R}^2}F_i\circ f_i\circ\Psi_i^*(x)\,{\rm d}x \nonumber\\ & \quad = \frac{1}{4\pi}\sum_{i=1}^l\int_{\mathbb{R}^2}\int_{\mathbb{R}^2}\ln\frac{k(x_{0,i})^\frac{1}{2}}{|x-y|}f_i\circ U_i(x)f_i\circ U_i(y)\,{\rm d}x\,{\rm d}y-\sum_{i=1}^l\int_{\mathbb{R}^2}F_i\circ f_i\circ U_i(x)\,{\rm d}x \nonumber\\ & \quad = \sum_{i=1}^l\frac{d_i^2}{8\pi}\ln k(x_{0,i})+C_0, \end{align}

where $C_0=({1}/{4\pi })\sum _{i=1}^l\int _{\mathbb {R}^2}\int _{\mathbb {R}^2}\ln ({1}/{|x-y|})f_i\circ U_i(x)f_i\circ U_i(y)\,{\rm d}x\,{\rm d}y-\sum _{i=1}^l\int _{\mathbb {R}^2}F_i\circ f_i\circ U_i(x)\,{\rm d}x$ is a constant dependent of $f_i$. Taking (4.36) into (4.35), we get the result.

Moreover, one can get the order of the functional $\mathcal {E}_i(\omega ^\varepsilon )$ and constants $\sigma _i^\varepsilon$ as follows

Lemma 4.13 For $\varepsilon$ sufficiently small, there holds

(4.37)\begin{align} \mathcal{E}_i(\omega^\varepsilon) & = \frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\frac{d_i^2}{8\pi}\ln k(x_{0,i})-\frac{d_i^2}{2}h(x_{0,i},x_{0,i})+C_1+o(1), \end{align}

(4.38)\begin{align} \sigma_i^{\varepsilon}& = \frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+ \frac{|d_i|}{4\pi}\ln k(x_{0,i})-|d_i|h(x_{0,i},x_{0,i})\nonumber\\ & \quad +sgn(d_i)\sum_{j=1,j\neq i}^ld_jG(x_{0,i},x_{0,j})+C_2+o(1), \end{align}

where $C_1=({1}/{4\pi }) \int _{\mathbb {R}^2}\int _{\mathbb {R}^2}\ln ({1}/{|x-y|})f_i\circ U_i(x)f_i\circ U_i(y)\,{\rm d}x\,{\rm d}y- \int _{\mathbb {R}^2}F_i\circ f_i\circ U_i(x)\,{\rm d}x$, $C_2=({1}/{2\pi })\int _{\mathbb {R}^2}\ln ({1}/{|x^*-y'|})f_i\circ U_i(y')\,{\rm d}y'$, and $x^*$ is any point of $\partial \{x\in \mathbb {R}^2\mid f_i\circ U_i(x)=0\}$.

Proof. Similarly as the proof of lemma 4.12, one can immediately get

\begin{align*} \mathcal{E}_i(\omega^\varepsilon) & =\frac{d_i^2}{4\pi}\ln\frac{1}{\varepsilon}+\frac{d_i^2}{8\pi}\ln k(x_{0,i})-\frac{d_i^2}{2}h(x_{0,i},x_{0,i})\\ & \quad +\frac{1}{4\pi} \iint_{\mathbb{R}^2\times \mathbb{R}^2}\ln\frac{1}{|x-y|}f_i\circ U_i(x)f_i\circ U_i(y)\,{\rm d}x\,{\rm d}y\\ & \quad - \int_{\mathbb{R}^2}F_i\circ f_i\circ U_i(x)\,{\rm d}x+o(1). \end{align*}

For (4.38), notice that for any $x\in \partial (supp(\omega _i^\varepsilon ))$, $\sigma _i^\varepsilon =sgn(d_i)\int _{D}G(x,\,y) \omega ^\varepsilon (y) \,{\rm d}\theta (y)$, which implies that

(4.39)\begin{align} \sigma_i^\varepsilon & =\frac{1}{2\pi}\int_{B_i}\ln\frac{1}{|x-y|}|\omega_i^\varepsilon(y)|\,{\rm d}\theta(y)-|d_i|h(x_{0,i},x_{0,i})\nonumber\\ & \quad +sgn(d_i)\sum_{j=1,j\neq i}^ld_jG(x_{0,i},x_{0,j})+o(1). \end{align}

Let $x=\varepsilon x'+\hat {X}^\varepsilon _i$, then $x'\in \partial (supp(\xi _i^\varepsilon ))$. So

(4.40)\begin{align} & \frac{1}{2\pi}\int_{B_i}\ln\frac{1}{|x-y|}|\omega_i^\varepsilon(y)|\,{\rm d}\theta(y) \nonumber\\ & \quad =\frac{|d_i|}{2\pi}\ln\frac{1}{\varepsilon}+\frac{1}{2\pi}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x'-y'|}k(\varepsilon y'+\hat{X}^\varepsilon_i)\xi_i^\varepsilon(y')\,{\rm d}y'. \end{align}

By lemma 4.10 and the continuity of $k$, we have

(4.41)\begin{align} & \frac{1}{2\pi}\int_{B_{\bar{ R}_1}(0)}\ln\frac{1}{|x'-y'|}k(\varepsilon y'+\hat{X}^\varepsilon_i)\xi_i^\varepsilon(y')\,{\rm d}y' \nonumber\\ & \quad =\frac{k(x_{0,i})}{2\pi}\int_{\mathbb{R}^2}\ln\frac{1}{|x'-y'|}f_i\circ\Psi_i^*(y')\,{\rm d}y'+o(1) \nonumber\\ & \quad =\frac{1}{2\pi}\int_{\mathbb{R}^2}\ln\frac{k(x_{0,i})^\frac{1}{2}}{|x^*-y'|}f_i\circ U_i(y')\,{\rm d}y'+o(1) \nonumber\\ & \quad =\frac{|d_i|}{4\pi}\ln k(x_{0,i})+\frac{1}{2\pi}\int_{\mathbb{R}^2}\ln\frac{1}{|x^*-y'|}f_i\circ U_i(y')\,{\rm d}y'+o(1), \end{align}

where $x^*$ is any point of $\partial \{x\in \mathbb {R}^2\mid f_i\circ U_i(x)=0\}$. Taking (4.40), (4.41) into (4.39), we get (4.38).

4.3. Proof of theorem 2.4

Proof. The proof is similar to the proof of theorem 2.3. By lemma 4.5, we know that $\omega ^\varepsilon$ has the form

\[ \omega^\varepsilon= \sum_{i=1}^lsgn(d_i)\frac{1}{ \varepsilon^2}f_i(\psi_i^{\varepsilon})\chi_{\{ \psi_i^{\varepsilon}>0 \}\cap B_i }, \]

where $\psi _i^{\varepsilon }=sgn(d_i)\mathcal {G}(k\omega ^{\varepsilon })-\sigma ^{\varepsilon }_i$ for some $\sigma ^{\varepsilon }_i\in \mathbb {R}.$ By lemma 4.6, $diam(supp(\omega _i^\varepsilon ))=O(\varepsilon )$.

Moreover, by proposition 4.7, the support set of $\omega _i^\varepsilon$ tends to $x_{0,i}$ as $\varepsilon \to 0^+$, namely

\[ \lim_{\varepsilon \to 0^+}\sup_{x\in supp(\omega_i^\varepsilon)}|x- x_{0,i}|=0. \]

By lemmas 4.10 and 4.13, we get (3)(4) in theorem 2.4. The proof of theorem 2.4 is thus complete.

4.4. Proof of theorems 1.3 and 2.2

Proof. Let $l=1$ and $d_1=I$ in theorem 2.4, we get theorem 2.2. Let $\lambda ={1}/{\varepsilon ^2}$ and $u^\lambda =u^\varepsilon,\, c^\lambda =-\mu ^\varepsilon$, we get theorem 1.3.