0$. We will write $p=p^*-\epsilon$. There are several motivations for studying this equation coming from geometry, analysis, mathematical biology, physics, etc. For $p=p^*$ and $\lambda=0$ equation \eqref{esfera} is related to the Yamabe problem and has been extensively studied (see for instance \cite{Struw}). It is well known that in this case the associated variational problem exhibits lack of compactness and that the Palais-Smale (PS) condition does not hold (\cite{LP},\cite{Struw}). Several authors have handled this difficulty by considering existence and blow up of solutions when $p\to p^*$ (\cite{LY}). In applications, equation \eqref{esfera} appears also in several contexts related to pattern formation. It is obtained when looking for steady state solutions of the reduced Gierer-Meinhardt system when the diffusion rate of the activator to that of the inhibitor tends to zero, which is known as the shadow equation \cite{NT}. Similar models appear in other contexts, such as flame propagation (see \cite{B}). Considering the problem on the sphere is also natural from several perspectives. Usually, pattern formation models are studied on plane domains with boundary; however, there are several organisms, for instance radiolarians, that exhibit (nearly) spherical symmetry. Some specimens of radiolarians show peaks similar to the ones depicted in figure 2. This was one of the main motivations for this work. In fact, it is perhaps more natural to visualize the dermis of an organism, on which the pattern formation process is taking place, as a surface, rather than a plane domain. It has been proposed that curvature effects can also be important in the selection of patterns and this makes it necessary to study the usual models on surfaces (\cite{PPS}). From the technical point of view, considering the sphere actually simplifies the problem, since boundary effects do not have to be taken into account. We can state our results in an informal way as follows (for a precise formulation see section~3): there are solutions which exhibit a prescribed number of peaks. Moreover, among solutions with a fixed number of concentrations (peaks), those with minimal energy are such that the centers of the concentrations tend to be as far as possible from each other, thus solving a sphere packing problem. A very similar result has been obtained by Gui and Wei \cite{GW} for domains in $\mathbb{R}^n$ with boundary. There are however several differences. First of all, these authors consider the equation \begin{gather*} \epsilon^2 \Delta u -u +f(u)=0,\quad u>0\quad \text{in } \Omega,\\ \frac{\partial u}{\partial \nu}=0 \quad \text{on } \partial\Omega \end{gather*} and their analysis is made taking $\epsilon$ as a small parameter. Moreover, as stated above, the presence of the boundary makes some of the estimates and technicalities somewhat involved. On the other hand, we take the coefficient $\epsilon^2$ of the Laplacian as fixed, and vary the exponent of the nonlinear term. Although these approaches are in some sense equivalent, as has been shown by Benci and Cerami (see for instance \cite{BC}), considering the exponent as our parameter allows us to use well known global compactness results (see \cite{Ba} or \cite{Struw}). The existence of $k$-peaked solutions (as defined in section 3) is obtained by minimization on a suitable space of functions, basically, those with the right $2\pi/k$-symmetry. We will construct an action that does not have fixed points on the sphere. Therefore, if a solution has a peak, by symmetry it has to have at least $k$. From the mathematical point of view, it is known that if $\lambda=0$ and $p=p^*$ there is a unique nontrivial positive solution (up to rotations and scalings), in other words, for the critical case, there are no $k$-peaked solutions on the sphere. Roughly speaking, this is due to the fact that the Palais-Smale condition is not valid any more and (PS) sequences that exhibit more than one peak are not compact. Using Bahri's terminology, these constitute critical points at infinity (\cite{Ba}). In our case, the fact that the exponent of the non linear term remains strictly smaller than the critical exponent implies that the (PS) condition still holds and therefore (PS) sequences with $k$ peaks do have subsequences converging to actual solutions. The analysis of existence and blow up when the linear term is not present has been the subject of much research; see for instance the work by Y. Li \cite{LY}, and references therein.\\ In order to make the reasoning more transparent, we include in section~3 an example in which the existence of multipeak solutions is more easily established. Namely, we look for solutions with two peaks. In this case, constraining the functional to even functions on the sphere is sufficient, since the action of the group, in this case the antipodal action, has only one fixed point, the origin, which does not belong to the sphere. The rest of the paper is organized as follows. In section~2 we formulate the variational problem and recall some well known facts. In section~3 we first state our results: existence of $k$-peaked solutions (theorem 3.1 and remark 3.1), and their geometric properties (theorem~3.2). In section 3.1 we give the proofs. We conclude in section~4 with some questions and open problems. \section{Formulation of the problem} We can work directly on $\mathbf{S}^n$ or use stereographic projection and work on $\mathbb{R}^n$, but generally computations are easier on $\mathbb{R}^n$ and many standard results are stated in that setting. Let $P=(0,\dots,0,1)$ be the north pole on $\mathbf{S}^n\subset \mathbb{R}^{n+1}$. Stereographic projection $\sigma:\mathbf{S}^n-\{P\} \to \mathbb{R}^n$ is defined by $\sigma(x^1,\dots,x^n,\xi)=(y^1,\dots,y^n)$ for $(x,\xi)\in \mathbf{S}^n-\{P\}$ where $$ y^j=x^j/(1-\xi). $$ We recall that equation \eqref{esfera} with $p=p^*-\epsilon$ is transformed, after stereographic projection on $\mathbb{R}^n$, into \begin{equation} -\Delta u+\frac{4\lambda}{(1+|y|^2)^2} u= \frac{(n-2)}{4(n-1)}\Big(\frac{2}{1+|y|^2}\Big)^{\frac{\epsilon}{2} (n-2)}u^p.\label{rn} \end{equation} Note that the term containing $d(n)$ cancels out. This equation is equivalent to \eqref{esfera} %the one on the sphere in the sense that if $v$ is a solution of \eqref{esfera} and $$ v(x)=\Big(\frac{2}{1+|y|^2}\Big)^{\frac{2-n}{2}} u(y) $$ (with $y=y(x)$) then $u$ solves equation (\ref{rn}) (see \cite{LP}). Associated with these two equations we have the following functionals $$ E^\epsilon_{\mathbf{S}^n}(v)=\int_{\mathbf{S}^n}\Big(\frac{1}{2}|\nabla v|^2 +\frac{(d(n)+ \lambda)}{2} v^2-\frac{v^{p^*+1-\epsilon}}{p^* +1-\epsilon}\Big)\,d\sigma $$ in $H^1(\mathbf{S}^n)$ and $$ E^\epsilon(u)=\int_{\mathbb{R}^n} \Big(\frac{1}{2}|\nabla u|^2+ \frac{2\lambda}{(1+|y|^2)^2}u^2- \frac{(n-2)}{4(n-1)} \big(\frac{2}{1+|y|^2}\big)^{\frac{\epsilon}{2} (n-2)} \frac{u^{p^*-\epsilon+1}}{p^*-\epsilon+1}\Big)\, dy $$ in $H^1(\mathbb{R}^n).$ Note that for small $\epsilon>0$ it is a standard fact that the Palais-Smale (PS) condition holds for $E^\epsilon_{\mathbf{S}^n}$ and for $E^\epsilon$ when restricted to bounded domains (\cite{Struw}) and solving equations \eqref{esfera} or \eqref{rn} is equivalent to finding critical points of the corresponding functionals in the following sense. The functional $E_{S^n}^{\epsilon}$ is not bounded below but one can add a constraint to be able to apply the direct method of the calculus of variations as done in [\cite{Struw}, I.2]. The idea is to consider the functional \begin{equation} \int_{\mathbf{S}^n}\big(\frac{1}{2}|\nabla v|^2+\frac{(d(n)+ \lambda)}{2} v^2\big)\,d\sigma \label{fcsn} \end{equation} restricted to the set where $\int_{S^n} |v|^{p^*+1-\epsilon}=1.$ Using Lagrange multipliers this yields the functional $S_{S^n}$ defined by: \begin{equation} S_{S^n}=\Big[ \int_{S^n} \frac 12 |\nabla v|^2 +\frac{(d(n)+\lambda)}{2} v^2 \Big] \Big/ \Big( \int_{S^n} |v|^{p^*+1-\epsilon}\Big) ^{2/(p^*+1-\epsilon)}.\label{2a} \end{equation} We will use a modification of a global compactness result by Struwe which we quote for completeness. For a bounded domain $\Omega\subset\mathbb{R}^n$ define $$ E_\lambda=\frac{1}{2}\int_\Omega \big(|\nabla u|^2+\lambda |u|^2\big)\, dx - \frac{1}{2^*}\int_\Omega |u|^{2^*}\,dx,$$ $2^*=2n/(n-2)=p^*+1$ and let $D^{1,2}(\Omega;\mathbb{R}^n)$ be the completion of $C_0^\infty(\Omega;\mathbb{R}^n)$ in the norm $\|u\|_{D^{1,2}}=\|\nabla u\|_{L^2}$. \begin{theorem}[{\cite[Thm. III. 3.1 p. 169]{Struw}}] \label{thm2.1} Suppose $\Omega$ is a bounded domain in $\mathbb{R}^n$, $n\ge 3$ and for $\lambda\in \mathbb{R}$ let $(u_m)$ be a Palais-Smale sequence for $E_\lambda$ in $H_0^{1,2}(\Omega)\subset D^{1,2}(\mathbb{R}^n)$. Then there exists an index $k\in {\bf N}_0$, sequences $(R^j_m), (x_m^j), \ 1\le j\le k$, of radii $(R^j_m)\to\infty(m\to\infty)$ and points $x^j_m\in \Omega$, a solution $u_0\in H_0^{1,2}(\Omega)\subset D^{1,2}(\mathbb{R}^n)$ to the problem \begin{equation} \label{1.1} \begin{gathered} -\Delta u =- \lambda u + u |u|^{2^*-2} \quad \text{in }\Omega\\ u>0 \quad \text{in } \Omega,\\ u=0 \quad\text{on } \partial\Omega \end{gathered} \end{equation} and non-trivial solutions $u^j \in D^{1,2}(\mathbb{R}^n)$, $1\le j\le k$, to the ``limiting problem'' $$ -\Delta u = u |u|^{2^*-2} \quad \text{in }\mathbb{R}^n, $$ such that a subsequence $(u_m)$ satisfies \begin{equation} \big\| u_m-u^0 - \sum_{j=1}^k u_m^j \big \|_{D^{1,2}(\mathbb{R}^n)} \to 0.\label{conv} \end{equation} Here $u_m^j$ denotes the rescaled function $$ u_m^j(x)=(R_m^j)^{{n-2\over 2}} u^j(R^j_m(x-x_m^j)), \quad 1\le j\le k,\; m\in {\bf N}. $$ Moreover $$ E_\lambda(u_m) \to E_\lambda(u^0)+\sum_{j=1}^k E_0(u^j). $$ \end{theorem} \begin{remark} \label{rmk2.1} \rm An analogous result when $\Omega$ is a compact subset of $\mathbf{S}^n$ (different from $\mathbf{S}^n$) can be found in \cite {Ba} or \cite{BaC}. We can in fact apply the result to $\mathbf{S}^n$ in this particular case because a minimizing sequence cannot blow-up at an infinite number of points on the sphere. Clearly any sequence with a finite number of blow-up points would have less energy. Therefore we can rotate so that the north pole is not a blow-up point. We can choose a ball $\tilde{ B}$ around the north pole where we can apply lemma 3.4 in \cite{MSW}: \end{remark} \begin{lemma}[{\cite[lemma 3.4]{MSW}}] \label{lem2.1} Let $\Omega_m=\{ x\in \mathbb{R}^n| \lambda_m^{-1/2} x\in \Omega\}$, $\Omega$ a bounded open set containing the origin, with $\lambda_m\to\infty$ as $m\to\infty$. Suppose $(v_m)_{m\in \mathbb{N}} \subset W^{1,2}(\Omega_m)$ is a sequence with $\|v_m\|_{W^{1,2}(\Omega_m)}$ uniformly bounded. If for any $R>0$, $$ \lim_{m\to \infty} \Big(\sup_{y\in\mathbb{R}^n} \int_{\Omega_m\cap B_R(y)} |v_m|^{p+1} \, dx\Big)=0, $$ then $$ \lim_{m\to\infty} \int_{\Omega_m} |v_m|^{p+1}\, dx=0. $$ \end{lemma} The proof of the above lemma can be found in \cite{W}. Since there is no blow-up in $\tilde{B}$, the condition of the lemma is satisfied and so $u_m\to 0$ strongly in $L^{p+1}(\Omega)$. By considering the complement $(\tilde B^\prime)^c$, where $\tilde B^\prime $ is a slightly smaller ball around the north pole, we can substitute the condition $u_m\in H_0^{1,2}(\Omega)$ in Theorem 2.1 by the fact that in a neighborhood of the boundary $u_m\to 0$ in the $L^{p+1}$ norm. \section{Main results} We begin by stating our results in a precise way. \begin{theorem} \label{thm3.1} Let $\lambda >0$ and $k\in\mathbb{N}^+$ be given. Then there exists an $\epsilon_0>0$ depending on $\lambda$ and $k$ such that equation \eqref{esfera} on $\mathbf{S}^n$, $n\ge 3$ and odd with $p=p^*-\epsilon$ has a positive solution, $u_\epsilon$, for all $0<\epsilon\le \epsilon_0$, which concentrates at $k$ different points $x^j\in \mathbf{S}^n,\ j=1,\dots k$. \end{theorem} \begin{remark} \label{rmk3.1} \rm That a positive solution, $u_\epsilon$, $0<\epsilon\le \epsilon_0$ concentrates at $k$ different points $x^j\in \mathbf{S}^n,\ j=1,\dots k$, or that $u_\epsilon$ has $k$ peaks means that when $\epsilon\to 0$, after rotations and rescalings, a subsequence converges strongly in $H^1(\mathbf{S}^n)$ to a linear combination of $k$ distinct solutions of the limiting problem \begin{equation} -\Delta_{\mathbf{S}^n} v + d(n) v = v^{p^*}.\label{lim} \end{equation} (Cf. Theorem~2.1, global compactness result.) \end{remark} Theorem 3.1 states that there is at least one solution with $k-$peaks. Now, for any fixed, small, positive $\epsilon$ the Palais-Smale condition is satisfied. Then, since solutions with exactly $k$ peaks and bounded energy form a compact set, we can assume that there is a solution $u^*_\epsilon$, (not necessarily the same one found in theorem 3.1) that, among all the $k$-peaked solutions, has least energy $S_{S^n}$. We now show that the behavior of the $k$ concentrations for $u^*_\epsilon$ when $\epsilon\to 0$ is governed by an underlying geometric problem: \begin{theorem} \label{thm3.2} Let $v_{\epsilon_i}$, with $\epsilon_i\to 0$ as $i\to\infty$ be a sequence of $k-$peaked solutions of equation \eqref{esfera} which have least energy $S_{S^n}$ as discussed in the preceding paragraph and let $x^j_i$, $j=1,\dots k,\ i=1,2,\dots$ be the center of the $j^{th}$ peak. Then there exists a subsequence (still denoted $x^j_i$) such that $x^j_i \to x^j_*$, where $x^j_*$ is a solution of the following maximization problem: $$ \max_{ x^j,x^l\in \mathbf{S}^n} \sum_{j\ne l} \| x^j-x^l\|^2 $$ where $x^j\ne x^l$ if $j\ne l$. \end{theorem} \begin{remark} \label{rmk3.2} \rm From theorem 3.2 it follows that the concentrations of least energy $k-$peaked solutions, as $\epsilon_i\to 0$, arrange themselves according to a packing problem on the sphere in the sense that the centers of the concentrations (peaks) tend to be as far away from each other as possible. \end{remark} \section*{Proofs} We begin by proving Theorem 3.1, the existence of $k-$peaked solutions. This is done by minimizing on a suitable set the functional $$ S_{S^n}(v)=\Big[ \int_{S^n} {1\over 2} |\nabla v|^2 +{(d(n)+\lambda) \over 2} v^2 \Big] \Big/\Big( \int_{S^n} |v|^{p^*+1-\epsilon}\Big) ^{2/(p^*+1-\epsilon)}. $$ We will restrict to functions which are invariant under the action of $g$ given by (\ref{g}) below. \subsection*{Case $k=2$} Before proving theorem 3.1 in the general case, we mention that the case $k=2$, namely existence of a solution with two peaks can easily be established. In fact, it is sufficient to minimize the corresponding functional, $S_{S^n}$, in the subspace of even functions, that is, functions invariant under the antipodal action $x\to -x$. Then by direct minimization and the fact that the (PS) condition holds for $1

0$ be given. Then there exists $t_0(\delta)$ and $C_1(n,k)$, $C_2(n,k)$, $e_1$, $e_2$ all greater than zero such that $$ I^\epsilon=C_1-C_2\sum_{j\ne l}^k \|y^j -y^l\| +e_1 +e_2 $$ where $e_1\to 0 $ when $\epsilon \to 0$, $e_2<\delta$ for all $t_0>t$. \end{proposition} \begin{proof} First observe that, in the limiting case, the quantity $$ J=\Big( \int_{\mathbf{S}^n} |\nabla \tilde{w}^j|^2\, d\sigma \Big)\Big/ \Big(\int_{\mathbf{S}^n} |\tilde{w}^j|^{p^*+1}\, d\sigma\Big)^{2/(p^*+1)} $$ does not depend on either $t$ or $y^j$. Indeed, \begin{equation} J=\frac{\int_{\mathbb{R}^n} \Big|\nabla \big(\frac{1/t}{1+|(\frac{z-z^j}{t})|^2} \big)^{(n-2)/2} \Big|^2\,|\mbox{Jac}|\, dz} {\Big(\int_{\mathbb{R}^n} \Big(\big( \frac{1/t}{1+\left(\frac{|z-z^j|}{t} \right)^2} \big)^{(n-2)/2}\Big)^{2n/(n-2)}\,|\mbox{Jac}|\, dz \Big)^{(n-2)/n} } \label{JR} \end{equation} and making the change of variables $ y=\frac{z-z^j}{t}$ expression \eqref{JR} becomes \begin{align*} &\int_{\mathbb{R}^n} |\nabla (\frac{1}{1+|y|^2})^{(n-2)/2}|^2 \frac{1}{t^2} t^n|\mbox{Jac}|\, dy \Big/ \Big( \int_{\mathbb{R}^n} \big(\frac{1}{1+|y|^2}\big)^n t^n |\mbox{Jac}|\, dy \Big)^{(n-2)/n}\\ &=\int_{\mathbb{R}^n} |\nabla (\frac{1}{1+|y|^2})^{(n-2)/2} | ^2|\mbox{Jac}|\, dy \Big/ \Big( \int_{\mathbb{R}^n} \big( \frac{1}{1+|y|^2}\big)^n \, |\mbox{Jac}|\, dy \Big)^{(n-2)/n}, \end{align*} which is independent of $t$ and $y^j$. By a similar change of variables we have that $$ \int_{\mathbf{S}^n } \frac{\lambda}{2} (\tilde{w}^j)^2\, d\sigma\Big/ \Big(\int_{\mathbf{S}^n} (\tilde{w}^j)^{p^*+1-\epsilon} \,d\sigma\Big)^{2/(p^* +1-\epsilon)} $$ tends to zero with $t$ uniformly in $0\le \epsilon\le 1$ because of the factor $t^{n-2n/(p^*+1-\epsilon)}$ which appears after scaling and has a positive exponent if $\epsilon<4/(n-2)$. We can take $t$ small enough so that the term containing $V_m^2$ in \eqref{I} is small, say less than $\epsilon_1$. Recall that the $L^p$ norm is an increasing and continuous function of $p$ and so, for $\epsilon $ small we can write for any $v\in L^{p^*+1}$, $$ \|v\|_{L^{p^*+1-\epsilon}}= \|v\|_{L^{p^*+1}}-p_1(\epsilon), $$ where $p_1(\epsilon)\ge 0$ and tends to zero with $\epsilon$. It also follows from Minkowski's inequality that $$ \Big( \int_{\mathbf{S}^n} v^{p^* +1}\,d\sigma \Big)^{1/(p^*+1)}= \sum_j\Big(\int_{\mathbf{S}^n} |\tilde{w}^j|^{p^*+1}\Big) ^{1/(p^*+1)} -p_2(t), $$ where $v(x)=v(\sigma^{-1}(y))=V(y)$ and $p_2(t)\ge 0$ and tends to 0 with $t$. So we can choose $t_0$ sufficiently small such that for all $\epsilon<\epsilon_0$ and $t