On the Neumann Problem of Hardy-Sobolev critical equations with the multiple singularities

Let $N \geq 3$ and $\Omega \subset \mathbb{R}^N$ be $C^2$ bounded domain. We study the existence of positive solution $u \in H^1(\Omega)$ of \begin{align*} \left\{ \begin{array}{l} -\Delta u + \lambda u = \frac{|u|^{2^*(s)-2}u}{|x-x_1|^s} + \frac{|u|^{2^*(s)-2}u}{|x-x_2|^s}\text{ in }\Omega\\ \frac{\partial u}{\partial \nu} = 0 \text{ on }\partial\Omega, \end{array}\right. \end{align*} where $0 < s <2$, $2^*(s) = \frac{2(N-s)}{N-2}$ and $x_1, x_2 \in \overline{\Omega}$ with $x_1 \neq x_2$. First, we show the existence of positive solutions to the equation provided the positive $\lambda$ is small enough. In case that one of the singularities locates on the boundary and the mean curvature of the boundary at this singularity is positive, the existence of positive solutions is always obtained for any $\lambda > 0$. Furthermore, we extend the existence theory of solutions to the equations for the case of the multiple singularities with different exponents.