Let \(n\ge 2\) , \(s\in (0,1)\) , and \(\Omega \subset \mathbb {R}^n\) be a bounded Lipschitz domain. This paper is devoted to the study of Green functions for the fractional Neumann problem \(\begin{aligned} \left\{ \begin{aligned} (-\Delta )^su&=f & \text {in}\ \ \Omega ,\\ \mathcal {N}_s u&=0 & \text {in}\ \ \mathbb {R}^n\setminus \Omega , \end{aligned}\right. \end{aligned}\) and the regional fractional Neumann problem \(\begin{aligned} \left\{ \begin{aligned} (-\Delta )^s_\Omega u&=f & \text {in}\ \ \Omega ,\\ \partial ^{2s-1}_{\varvec{\nu }} u&=0 & \text {on}\ \ \partial \Omega , \end{aligned}\right. \end{aligned}\) often referred to as Neumann functions. Here, \(\mathcal {N}_s\) and \(\partial ^{2s-1}_{\varvec{\nu }}\) denote the nonlocal normal derivative in \(\mathbb {R}^n\setminus \Omega \) and the fractional normal derivative on \(\partial \Omega \) respectively. Precisely, we establish the existence and uniqueness of Neumann functions for the fractional Neumann problem and the regional fractional Neumann problem respectively. Moreover, we also obtain global pointwise upper bound estimates, fractional Sobolev type estimates, and Hölder continuity estimates for such Neumann functions. Furthermore, a representation formula for the weak solution to the nonhomogeneous (regional) fractional Neumann problem is obtained by using the Neumann function.
