  Let $\Omega$ be a bounded domain in $\mathbb{R}^2$ with smooth boundary, we
study the following elliptic Dirichlet problem $$ \begin{cases}
-\Delta\upsilon= e^{\upsilon}-s\phi_1-4\pi\alpha\delta_p-h(x)\,\,\,\,
\,\textrm{in}\,\,\,\,\,\Omega,\\[2mm] \upsilon=0 \quad\quad\quad\quad\quad\quad
\qquad\qquad\quad\quad\,\,\,\, \textrm{on}\,\ \,\partial\Omega, \end{cases} $$
where $s>0$ is a large parameter, $h\in C^{0,\gamma}(\overline{\Omega})$,
$p\in\Omega$, $\alpha\in(-1,+\infty)\setminus\mathbb{N}$, $\delta_p$ denotes
the Dirac measure supported at point $p$ and $\phi_1$ is a positive first
eigenfunction of the problem $-\Delta\phi=\lambda\phi$ under Dirichlet boundary
condition in $\Omega$. If $p$ is a strict local maximum point of $\phi_1$, we
show that such a problem has a family of solutions $\upsilon_s$ with arbitrary
$m$ bubbles accumulating to $p$, and the quantity
$\int_{\Omega}e^{\upsilon_s}\rightarrow8\pi(m+1+\alpha)\phi_1(p)$ as
$s\rightarrow+\infty$.
