Spiral neighbourhoods
R, 156 bytes



function(n){g=function(h)c(0,cumsum(h((4*(0:(n+2)^2)+1)^.5%%4%/%1/2)))
x=g(sinpi)
y=g(cospi)
a=x[n]
b=y[n]
which(x==a&amp;(y==b+1|y==b-1)|y==b&amp;(x==a+1|x==a-1))}


Try it online!


posted another R answer since it's a slightly different approach than @ngn
1-indexed
neighbours are always sorted by ascending value
saved 6 bytes removing round and using cospi(x)/sinpi(x) which are more precise than cos(x*pi)/sin(x*pi) in case of half numbers (0.5, 1.5 etc...)
saved another byte removing the minus on y coordinates since the result is the same (just up/down neighbours are reversed)




Explanation :

If we look at the matrix coordinates of the values, considering the first value 0 placed at x=0, y=0, they are :

x = [0,  1,  1,  0, -1, -1, -1,  0,  1,  2,  2,  2,  2,  1,  0, ...] 
y = [0,  0,  1,  1,  1,  0, -1, -1, -1, -1,  0,  1,  2,  2,  2, ...]


The x coordinates follow the A174344 OEIS sequence with the recursive formula :  

a(1) = 0, a(n) = a(n-1) + sin(mod(floor(sqrt(4*(n-2)+1)),4)*pi/2)


The same formula holds for y matrix coordinates, but with cos instead of sin and negated :

a(1) = 0, a(n) = a(n-1) - cos(mod(floor(sqrt(4*(n-2)+1)),4)*pi/2)


So, in R we can translate the formula to this function, taking sinpi/cospi as parameter :

g=function(h)c(0,cumsum(h((4*(0:(n+2)^2)+1)^.5%%4%/%1/2)))


and we generate the two coordinates vectors (we don't negate the y coords since we'll get the same result, just with up/down neighbours reversed) :

x=g(sinpi)
y=g(cospi)


Note that we have generated (n+2)^2 coordinates, which are more than the minimum necessary coordinates containing both n and their neighbours (a tighter bound would be (floor(sqrt(n))+2)^2 but unfortunately is less "golfy").

Therefore, now that we have all the coordinates, we first search the coordinates a,b corresponding to our n : 

a=x[n]
b=y[n]


finally we select the positions of their neighbours, i.e. : 


the up/down neighbours where x == a and y == b+1 or b-1
the right/left neighbours where y == b and x == a+1 or a-1


using :

which(x==a&amp;(y==b+1|y==b-1)|y==b&amp;(x==a+1|x==a-1))