Let $f(x) \in \mathbb{Z}[x]$ be a degree $d>1$ polynomial with integer coefficients. Define

$$r(n,M) := | \{x,y \in \mathbb{Z} : f(x)+f(y) = n \}|. $$

My question is:

Does the estimate $r(n,M) \ll_{\epsilon} n^{\epsilon} $ hold for every $\epsilon >0$?

In certain cases (such as $ f(x)=x^{2k}$), one can "factor" the problem and deduce the desired result from the divisor bound. However, I do not see how to approach the general case in this manner. I am aware that there is a weaker but more general result of Bombieri and Pila which states that

$$r'(n,M) := | \{x,y \in \mathbb{Z} : f(x,y) = n, |x|,|y| < M \}| $$

satisfies $r'(n,M) \ll_{\epsilon} M^{1/d + \epsilon} $when $f(x,y)$ is an absolutely irreducible polynomial of degree $d$. In this greater level of generality this is nearly best possible as can be seen by taking $f(x,y) = x^d -y$.