Let \(s\in (0,1)\) , \(\varepsilon >0\) and let \(\Omega \) be a bounded smooth domain. Given the problem \(\begin{aligned} \varepsilon ^{2s}(-\Delta )^{s} u + V(x)u = |u|^{p-1}u \quad \text{ in } \; \Omega , \end{aligned}\) with Dirichlet boundary conditions and \(1<p<(n+2s)/(n-2s)\) , we analyze the existence of positive multi-peak solutions concentrating, as \(\varepsilon \rightarrow 0\) , to one or several points of \(\Omega \) . Under suitable conditions on V, we construct positive solutions concentrating at any prescribed set of its non degenerate critical points. Furthermore, we prove existence and non existence of clustering phenomena around local maxima and minima of V, respectively. The proofs rely on a Lyapunov-Schmidt reduction where three effects need to be controlled: the potential, the boundary and the interaction among peaks. The slow decay of the associated ground-state demands very precise asymptotic expansions.
