A polynomial formula for the primes (with 26 variables) was presented by Jones, J., Sato, D., Wada, H. and Wiens, D. (1976). Diophantine representation of the set of prime numbers. American Mathematical Monthly, 83, 449-464.
The set of prime numbers is identical with the set of positive values taken on by the polynomial
$(k+2)(1-(wz+h+j-q)^2-((gk+2g+k+1)\cdot(h+j)+h-z)^2-(2n+p+q+z-e)^2-(16(k+1)^3\cdot(k+2)\cdot(n+1)^2+1-f^2)^2-(e^3\cdot(e+2)(a+1)^2+1-o^2)^2-((a^2-1)y^2+1-x^2)^2-(16r^2y^4(a^2-1)+1-u^2)^2-(((a+u^2(u^2-a))^2-1)\cdot(n+4dy)^2+1-(x+cu)^2)^2-(n+l+v-y)^2-((a^2-1)l^2+1-m^2)^2-(ai+k+1-l-i)^2-(p+l(a-n-1)+b(2an+2a-n^2-2n-2)-m)^2-(q+y(a-p-1)+s(2ap+2a-p^2-2p-2)-x)^2-(z+pl(a-p)+t(2ap-p^2-1)-pm)^2)$
as the variables range over the nonnegative integers.
I am asking if there exist a similar result for the primes of the form $n²+1$, i.e., this set of prime numbers is identical with the set of positive values taken on by certain polynomial.