We study the following Neumann boundary problem associated with the stationary solutions of the Keller-Segel system, a basic model of chemotaxis phenomena: \(\begin{aligned} \left\{ \begin{aligned} -\Delta _g u +\beta u&=\lambda \left( \frac{Ve^u}{\int _{\Sigma } Ve^u dv_g}-1\right) & \text{ in } \mathring{\Sigma }\\ \partial _{ \nu _g} u&=0 & \text{ on } \partial \Sigma \end{aligned}\right. , \end{aligned}\) on a compact Riemann surface \((\Sigma , g)\) of unit area, with interior \(\mathring{\Sigma }\) and smooth boundary \(\partial \Sigma \) . Here, \(\Delta _g\) denotes the Laplace-Beltrami operator, \(dv_g\) is the area element of \((\Sigma , g)\) , \(\nu _g\) is the unit outward normal to \(\partial \Sigma \) , \(\lambda \) and \(\beta \) are non-negative parameters, and V is non-negative with finite zero set. For any \(m\in \mathbb {N}\) and \(k,l\in \mathbb {N}\cup \{ 0\}\) with \(m=2k+l\) , we establish a sufficient condition on V for the existence of a sequence of blow-up solutions concentrating at k points in the interior and l points on the boundary as \(\lambda \) approaches the resonant values \(4\pi m\) . Moreover, our analysis extends to the corresponding singular problem.
