Find all positive integers $n$ ,such that the polynomial $ a_1^n+a_2^n + ... + a_n^n - n a_1 a_2 ... a_n $is irreducible over the real numbers

Short answer : Your conjecture is true, and the polynomial is in fact irreducible over $\mathbb C$. This is because your polynomial is basically a cyclotomic polynomial plus a negligible pertubation. Your polynomial is so close to a cyclotomic one that if it had any nontrivial divisor, it would have to be a "cyclotomic" divisor, but the perturbation prevents that.

Detailed answer :

Let $P(a_1,a_2,\ldots,a_n)=\sum_{k=1}^n a_k^n-na_1a_2\ldots a_n$ be your polynomial, and let $P^{\tilde{}}(X,Y)=P(1,1,\ldots,1,X,Y)$. It will suffice to show that $P^{\tilde{}}$ is irreducible in ${\mathbb C}[X,Y]$ (indeed, suppose that $P=A\times B$ is a nontrivial factorization of $P$. Then $P^{\tilde{}}=A^{\tilde{}} \times B^{\tilde{}}$, so one of the two factors (say $A^{\tilde{}}$) is constant in ${\mathbb C}[X,Y]$. Then $B^{\tilde{}}$ has maximal degree $n$ in both $X$ and $Y$, which implies the same for $B$. So $A$ has degree $0$ in both $X$ and $Y$. Then $A$ must divide all the $(X,Y)$-coefficients in $P$, namely $\sum_{k=1}^{n-2} a_k^n, -na_1a_2\ldots a_{n-2},1$ and $1$. So $A$ must be constant).

By some asymptotic root analysis (see appendix), there are constants $M,c>0$ such that for any $x\in [M,\infty)$, the polynomial $P^{\tilde{}}(x,Y)$ has roots $\alpha_1(x),\ldots,\alpha_n(x)$ with $$|\alpha_k(x)-\zeta_k x|\leq\frac{c}{x^{n-3}}\tag{1}$$ for each $k$, where $\zeta_k=e^{\frac{k}{n}i\pi}$.

Denote by $s_k$ the $k$-th elementary symmetric polynomial, $s_k(u_1,\ldots,u_n)=\sum_{i_1<\ldots i_k}u_{i_1}\ldots u_{i_k}$. And for $I\subseteq \lbrace 1,2,\ldots, n\rbrace$, say $I=\lbrace i_1,i_2,\ldots i_r \rbrace$ put $s_I(u_1,\ldots,u_n)=s_{r}(u_{i_1},u_{i_2},\ldots,u_{i_r})$. Finally let $\alpha(x)=(\alpha_1(x),\ldots,\alpha_n(x))$ and $\zeta=(\zeta_1,\ldots,\zeta_n)$. By multiplying and adding, we have for any $I\subseteq \lbrace 1,2,\ldots, n\rbrace$, a constant $c_I>0$ such that $$ |s_I(\alpha(x))-s_I(\zeta)x^{|I|}|\leq \frac{c_I}{x^{n-(|I|+2)}} \tag{2} $$

Suppose that there is a nontrivial factorization $P^{\tilde{}}=A \times B$ in ${\mathbb C}[X,Y]$. We may suppose that $A$ and $B$ are monic in $Y$. Note that of the two factors (say $A$) has degree $\leq n-3$ in $Y$, otherwise both factors have degree at least $n-2$, whence $n\geq 2(n-2)$ which is false for $n\geq 5$ (the case $n=4$ is similar and simpler).

Let $r$ be the degree of $A$ in $Y$. For any $x$, $A(x,Y)=\prod_{i\in I_x}(Y-\alpha_i(x))$ for some $I_x\subseteq \lbrace 1,2,\ldots, n\rbrace$ with cardinality $r$. Since there are only finitely many possible values for $I_x$, by the pigeon-hole principle there is an infinite $X\subseteq [M,\infty)$ and a fixed $I$ such that $I_x=I$ for all $x\in X$. Then for any $J\subseteq I$, $s_J(\alpha(x))$ is a coefficient appearing in $A(x,Y)$, so it must be a polynomial in $x$. Then (2) (used with $J$ in place of $I$) forces this polynomial to be $s_J(\zeta)x^{|J|}$. It follows that for any $x\in X$, $A(x,Y)= \prod_{i\in I}(Y-\zeta_ix)$. This identity holds for any $x$ because $X$ is infinite. So $A$ divides $X^n+Y^n$. But then $A$ also divides $P^{\tilde{}}-(X^n+Y^n)=-nXY+(n-2)$, which forces $r=1, A=Y-\zeta_kX$ for some $k$, and is clearly impossible. This finishes the proof.

Appendix : asymptotic root analysis

Here when we say "constant" we mean "independent of $x$". Let

$$ Q_x(T)=P^{\tilde{}}(x,xT)=x^n(T^n+1)-nx^2T+n-2 $$

Step 1 : There are two constants $K_1,M_1>0$ such that $|t| \leq K_1$ whenever $t$ is a root of $Q_x$ and $x\in (M_1,\infty)$.

It suffices to take $K_1,M_1$ such that $M_1\geq 1,K_1^{n-1}\geq n, K_1^n\geq 6(n-2)$. Indeed, if $|t|>K_1$ and $x\in (M_1,\infty)$, we have $|x^nt^{n-1}| \geq x^n K_1^{n-1}$ whence $|x^nt^{n-1}-nx^2| \geq x^n K_1^{n-1} -nx^2 \geq x^n\frac{K_1^{n-1}}{2}$ and $|x^nt^n-nx^2t|\geq x^n\frac{K_1^n}{2}$, and eventually $|Q_x(t)| \geq x^n\frac{K_1^n}{2}-((n-2)+x^n) \geq x^n\frac{K_1^n}{6}$, so $t$ cannot be a root of $Q_x$.

Step 2 : There are two constants $K_2,M_2>0$ such that $|t^n+1| \leq \frac{K_1}{x^{n-2}}$ whenever $t$ is a root of $Q_x$ and $x\in (M_2,\infty)$.

Thanks to step 1 and the fact that $t^n+1=\frac{nx^2t+n-2}{x^n}$, we see that we can take $K_2=n(K_1+1),M_2=M_1$.

The inequality in step 2 can be rewritten $\prod_{k=1}^n |t-\zeta_k| \leq \varepsilon$ where $\varepsilon =\frac{K_1}{x^{n-2}}$. In particular, $t$ satisfies at least one of the $n$ inequalities $|t-\zeta_k| \leq \varepsilon^{\frac{1}{n}}$. For large enough $x$, those inequalities will be mutually exclusive, showing

Step 3 : There are two constants $K_3,M_3>0$ such that whenever $x\in (M_3,\infty)$, $Q_x$ has roots $t_1,\ldots,t_n$ such that $|t_k-\zeta_k| \leq \frac{K_3}{x^{\frac{n-2}{n}}}$.

Putting $\varepsilon=t_k-\zeta_k$, we then have $\varepsilon \to 0$ as $x\to \infty$. Let $C=\frac{(\zeta+\varepsilon)^n-\zeta^n-n\zeta^{n-1}\varepsilon}{\varepsilon^2}$. By expanding the Newton binomial, we have that $C \to \frac{n(n-1)}{2}\zeta^{n-2}$ when $x\to \infty$. Now, note that $0=Q_x(t)=Q_x(\zeta_k+\varepsilon)$ can be rewritten as $$ \varepsilon=\frac{nx^2(\zeta+\varepsilon)-(n-2)+C\varepsilon^2}{n\zeta_k^{n-1}x^n} $$

We deduce

Step 4 : There are two constants $K_4,M_4>0$ such that whenever $x\in (M_4,\infty)$, $Q_x$ has roots $t_1,\ldots,t_n$ such that $|t_k-\zeta_k| \leq \frac{K_4}{x^{n-2}}$.

Multiplying by $x$ in step 4, we obtain (1) as announced.