Boundary concentration phenomena for an anisotropic Neumann problem in $\mathbb{R}^2$

Given a smooth bounded domain $\Omega$ in $\mathbb{R}^2$, we study the following anisotropic Neumann problem $$ \begin{cases} -\nabla(a(x)\nabla u)+a(x)u=\lambda a(x) u^{p-1}e^{u^p},\,\,\,\, u>0\,\,\,\,\, \textrm{in}\,\,\,\,\, \Omega,\\[2mm] \frac{\partial u}{\partial\nu}=0\,\, \qquad\quad\qquad\qquad\qquad\qquad\qquad \ \ \ \ \,\qquad\quad\, \textrm{on}\,\,\, \partial\Omega, \end{cases} $$ where $\lambda>0$ is a small parameter, $0<p<2$, $a(x)$ is a positive smooth function over $\overline{\Omega}$ and $\nu$ denotes the outer unit normal vector to $\partial\Omega$. Under suitable assumptions on anisotropic coefficient $a(x)$, we construct solutions of this problem with arbitrarily many mixed interior and boundary bubbles which concentrate at totally different strict local maximum or minimal boundary points of $a(x)$ restricted to $\partial\Omega$, or accumulate to the same strict local maximum boundary point of $a(x)$ over $\overline{\Omega}$ as $\lambda\rightarrow0$.