Question: For any $a\in \Bbb R$ and $n\in\Bbb N^+$, prove all the roots of $$(x+a)^n+x^n+(x-a)^n=0$$ are either real or purely imaginary.
Some thoughts so far:
If $z$ is a root of the equation, it is easy to verify that $-z,\overline z, -\overline z$ are the roots, too. I tried to put $z=u+vi$ and $\overline z=u-vi$ into the equaion but it looks more complicated.
Any ideas to proceed?
Your conjecture is true. In fact, we will prove the following:
Theorem: Let $P_{n,a}(x)$ be the polynomial $$P_{n,a}(x) = (x-a)^n + x^n + (x+a)^n,$$ where $n\in \mathbb{N^+}$ and $0\neq a\in \mathbb{R}$. Then all roots of $P_{n,a}$ lie on the imaginary axis.
Proof: We make a few simplifications first. First note that the theorem holds for $P_{n,a}$ if and only it holds for $P_{n,1}$ since $$(x-a)^n + x^n + (x+a)^n = 0 \iff \left(\frac{x}{a} - 1\right)^n + \left(\frac{x}{a}\right)^n + \left(\frac{x}{a} + 1\right)^n = 0,$$ so we will work with $P_n \equiv P_{n,1}$ without loss of generality from now on.
Now, note that $x=0$ is a root if and only if $n$ is odd. Because of this, it will also make things slightly more convenient for us to assume that $n$ is even. We can also do this without loss of generality because of the Gauss-Lucas Theorem and the fact that $P'_n(x) = nP_{n-1}(x)$. So it suffices to work with $P_{2n}$.
Finally, with even exponent it suffices to consider the equivalent equation $$x^{2n}\left(1-\frac{1}{x}\right)^{2n} + x^{2n} + x^{2n}\left(1+\frac{1}{x}\right)^{2n} = 0 \iff (1-iy)^{2n} + 1 + (1+iy)^{2n} = 0,\tag{*}$$ where we write $iy = x^{-1}$.
Let $z = 1+iy$. The solutions to $(*)$ with $y$ purely real (or equivalently, $x$ purely imaginary) corresponds to solutions of the equation $$\mathrm{Re}(z^{2n}) = -\frac{1}{2}.\tag{**}$$ Now, assume $y$ is real. We will show that $(**)$ has $2n$ solutions, and hence that $P_{2n}(x) = 0$ has $2n$ imaginary roots.
Let us write $z = re^{i\theta}$ where $r=\sqrt{1+y^2}$ and $\tan \theta = y$. Then $$z^{2n} = (1+y^2)^ne^{2ni\theta},$$ or equivalently $$\mathrm{Re}(z^{2n}) = (1+\tan^2 \theta)^n\cos(2n\theta) = \frac{\cos(2n\theta)}{\cos^{2n}(\theta)}.$$ We can write $\cos(2n\theta)$ in terms of the Chebyshev polynomials and hence we are looking for solutions of $$T_{2n}(\cos\theta) + \frac{1}{2}\cos^{2n}\theta = 0,$$ therefore we will be finished if we can show that $$T_{2n}(x) + \frac{1}{2}x^{2n} = 0 \tag{***}$$ has $n$ distinct solutions in the interval $(0,1)$, which will correspond to $2n$ distinct values of $\theta \in (-\pi/2 , \pi/2)$ which parametrizes the solutions of $(**)$.
Now, one of the defining features of Chebyshev polynomials is that all of their roots lie within the interval $(-1,1)$ and are symmetric about $0$. Hence $T_{2n}$ has precisely $n$ distinct positive real roots in $(0,1)$. Moreover, the magnitudes of the local maxima/minima of $T_n$ within the interval $(-1,1)$ are precisely $\pm 1$.
The addition of $\frac{1}{2}x^{2n}$ corresponds to a perturbation to $T_{2n}$ of magnitude less than $1$ on $(0,1)$ and hence cannot change the signs of any of local maxima/minima. Therefore by the intermediate value theorem, the addition of $\frac{1}{2}x^{2n}$ cannot change the fact that there exists $n$ distinct roots in $(0,1)$. This proves the result. $\square$