We are concerned with the study of \(C^{2, \gamma }_{loc}\) -positive solutions to the Choquard equation \( \displaystyle -\Delta u=\lambda (|x|^{-\alpha }*u^p)u^q \quad \text{ in } \mathbb {R}^N_+=\mathbb {R}^{N-1}\times (0, \infty )\, , N\ge 2, \) subject to the Robin boundary condition \(\beta u+\frac{\partial u}{\partial \nu }=g\) on \(\partial \mathbb {R}^N_+\) . Here \(p,q, \lambda >0\) , \(\alpha \in (0,N)\) , \(\beta \in \mathbb {R}\) and \(0\le g\in C^{1, \gamma }_{loc}(\mathbb {R}^{N-1})\) , \(\gamma \in (0,1)\) . First we show that if \(\beta <0\) then no solutions exist. Further, in the case \(\beta =0\) we prove that the existence of a \(C^{2, \gamma }_{loc}\) solution is closely related to the rate at which g decays at infinity. To this aim, we assume \(g(x')\simeq (1+|x'|)^{-m}\) , \(m>0\) , and determine the exact range of \(p, q, \alpha , m\) for which a \(C^{2, \gamma }_{loc}\) solution exists. A similar discussion arises in the case \(\beta >0\) . Our approach combines integral estimates and integral representation of superharmonic functions together with a new sub and supersolution method devised in a nonlocal setting.
