Is the "cyclotomic diagonalization" always squarefree? For every integer $n\ge 2$ , define $$f(n):=\Phi_n(n)$$ where $\Phi(n)$ is the $n$ th cyclotomic polynomial. This can be considered as the "cyclotomic diagonalization"

Prove or disprove the conjecture that $f(n)$ is squarefree for every integer $n\ge 2$

The partial factorizations I determined with PARI/GP (usually, the factors upto $12-13$ digits are found) gave no exponent greater than $1$ upto $n=200$. A rigorous test upto $n=46$ revealed that upto this limit there are only squarefree values.
If it cannot be proven, a good heuristic would also be welcome.
 A: We want to show that the problem is equivalent to
For any number $n$ and any prime $q \mid n^n-1$, then $ord_{q^2}(n) \neq n$
We will see at the end why we expect this to be false, but also that counterexamples are very rare and are likely to have more than forty digits.
Recall the Moebius function defined here, and the notation $v_p(k)$ for the maximum exponent $\alpha$ such that $p^{\alpha} \mid k$. Notice that $v_p(xy) = v_p(x)+v_p(y)$ behaves like a logarithm. The key observation is the following:
Lemma. It holds
$$ (**) \ \ \ \ v_p(\Phi_n(n) ) = \sum_{d \mid n} \mu \left (\frac{n}{d} \right ) v_p(n^d-1) $$
Proof. It is well known that
$$ \prod_{ d \mid n} \Phi_d(x) = x^n -1 $$
Now for any natural number $x$, taking $v_p$ yields
$$ \sum_{d \mid n} v_p \Phi_d(x) = v_p(x^n-1) $$
Applying the Moebius Inversion Formula we get
$$ v_p \Phi_n(x) = \sum_{d \mid n} \mu \left ( \frac{n}{d} \right ) v_p(x^d-1)$$
and the thesis follow by specializing at $x=n$.
Main Proof. Now that we have a formula for the exponent, let's compute it. Set $r:= ord_p(n)$. Note that if $r$ does not divide $n$, it does not divide any of its divisors $d \mid n$; it follows that $n^d \neq 1 \pmod{p}$ for all $d$. Thus equation (**) yields $v_p(\Phi_n(n) ) = 0$. Let us suppose for the rest that $r \mid n$ and set $n'=n/r$.
Case 1. $n'=1$, that is $ord_p(n) = n$. In this case the only summand that contributes to the summation in (**) is $v_p(n^n-1)$. This is $\le 1$ if and only if $p^2 \nmid n^n - 1$. Note that $n = ord_p(n) \mid ord_{p^2}(n) $:

*

*If the $ord_{p^2}(n)$ is equal to $n$, then $p^2 \mid n^n -1 $ and we found a counterexample;

*If instead $ord_{p^2}(n)>n$, we must have that $n^n \neq 1 \pmod{p^2}$ and we are happy.

Case 2 $n' >1$. This part can seem heavy, but if you are familiar with Moebius Function you will recognize that this is just ordinary Yoga.
The summands that contribute to (**) are all the $d$ that are divisible by $r$. Note that any number $d$ such that $r \mid d \mid n$ can be written as $rd'$, where $d' \mid n'=n/r $. With this observation in mind we write
$$ v_p(\Phi_n(n) ) = \sum_{r \mid d \mid n} \mu \left (\frac{n}{d} \right ) v_p(n^d-1) = 
\sum_{d' \mid n'} \mu \left (\frac{r n'}{r d'} \right ) v_p(n^{rd'}-1) = \sum_{d' \mid n'} \mu \left (\frac{n'}{d'} \right ) v_p((n^r)^{d'}-1)
$$
Now we use the Lifting The Exponent Lemma to split the $v_p$ factor. We get
$$ \sum_{d' \mid n'} \mu \left (\frac{n'}{d'} \right )( v_p(n^r-1) + v_p(d') )= v_p(n^r-1) \sum_{d' \mid n'} \mu \left (\frac{n'}{d'} \right ) + \sum_{d' \mid n'} \mu \left (\frac{n'}{d'} \right ) v_p(d') $$
Since $n'>1$, the first summand is zero: indeed $\sum_{k \mid h} \mu(h/k) = 0$ for $h> 1$. The second factor is more interesting. Suppose $n'= p^{\alpha} m$ with $(m,p) =1$. Every divisor of $n'$ will uniquely write as $p^{\beta} m'$, where $\beta \le \alpha$ and $m' \mid m$. Using the fact that $\mu(ab) = \mu(a)\mu(b)$ for $a,b$ coprimes we get
$$\sum_{d' \mid n'} \mu \left (\frac{n'}{d'} \right ) v_p(d') = \sum_{\beta=0}^{\alpha} \sum_{m' \mid m} \mu(p^{\alpha-\beta}) \mu(m/m') \beta = \sum_{\beta= 0}^{\alpha} \mu(p^{\alpha-\beta})\beta \sum_{m' \mid m} \mu(m/m')$$
We have to do a little distinction here. If $m>1$ this is zero again by the "summation moebius" property. If $m=1$, we are basically saying that $n'=p^{\alpha}$. In this case the formula reduce to
$$ \sum_{\beta=0}^{\alpha} \mu(p^{\alpha-\beta}) \beta = \beta - (\beta-1) = 1$$
And... Voilà!
Heuristics. We want to show that it is almost certain that a counterexample exist, but it is very likely to have hundreds of digits. Define for $d \mid p(p-1)$:
$$ X(p) = \{ 1 \le z \le p^2: \ \ (z,p) =1 \} $$
$$ X_d(p) := \{1 \le z \le p^2 : \ \ (z,p) = 1 \text{ and } ord_{p^2}(z) = d \} $$
Note that
$$\# X(p) = \phi(p^2) = p(p-1) \ \ \ \ \#X_d(p) = \#\{\text{roots of } x^d-1 \text{ in } \mathbb{F}_p\} = d$$
It is reasonable to assume that $X_d(p)$ is a random set in $X(p)$; indeed, if I am not wrong, there are some results of Davenport, Erdos, Vinogradov showing that quadratic residues, cubic residues ecc. are equidistributed among $X(p)$, and you can say that $ord_{p^2}(z) = k$ if and only if $z$ is a $\phi(p)/k$ residue but not a $\phi(p)/kq$ residue for all $q \mid \phi(p)/k$. However, I am not expert enough to produce a real proof.
We say that $p$ is a counterexample to the cdp (cyclotomic diagonalization problem) if $p^2$ appears as a factor of some $\Phi_n(n)$. We have shown above that this can happen iff $ord_{p^2}(n) = n$. Note that such an $n$ should be a divisor of $\phi(p) = p(p-1)$, because $n = ord_{p^2}(n) \mid p(p-1)$ and must be also coprime with $p$ (otherwise $p \nmid n^n-1$). On balance, the $n$ that make $p$ a counterexample to the cdp must be searched among the divisors of $(p-1)$.
For any such divisor $b$, we must hope that $b \in X_b(p)$. Since $X_b(p)$ is a random set in $X(p)$, this happens with probability $\frac{b}{p(p-1)}$. So the probability that this does not happen for any divisor $b$ is
$$ \prod_{b \mid (p-1) } \left ( 1- \frac{b}{p(p-1)} \right ) $$
The logarithm of this sum can be estimated as
$$ - \sum_{b \mid (p-1) } \frac{b}{p(p-1)} = - \frac{1}{p(p-1)} \sum_{b \mid q-1} b \le - \frac{1}{p} $$
Using the fact that $\sum_{b \mid q-1} b \ge q-1$. So the probability that $p$ is a counterexample is $\ge $
$$ 1- e^{-1/p} \simeq \frac{1}{p} $$
We use this logarithm-product strategy to estimate the probability $Q_x$ of finding a counterexample $\le x$. The probability of not having counterexamples up to $x$ is
$$ 1- Q_x = \prod_{p \le x} \left ( 1 - \frac{1}{p} \right ) \simeq exp(-\sum_{p \le x} \frac{1}{p} )\simeq exp(-\log \log x) = \frac{1}{\log(x)}$$
From which we get $Q_x \simeq 1- \frac{1}{\log x} $. This tends to $1$, so it is almost certain that a counterexample exists. On the other side, it goes to zero very slowly.
Example. Suppose you want to be $99\%$ sure to find a counterexample up to $x$. Then
$$\frac{99}{100} \simeq 1- \frac{1}{\log x} \ \ \Leftrightarrow \ \ \ \log x \simeq 100 $$
Which implies $x \simeq e^{100} \simeq 10^{43}$.
P.S. I am missing the case $p=2$, in which the Lifting the Exponent Lemma is slightly different; it should work the same way, though (the only different passage is the last, I which I think it is enough to verify the case $\Phi_{2^n}(x) = x^{2^{n-1}}+1$ :) )
A: Checking the pairs $n,p\le 10^6$ , I found one counterexample , namely :
$$283411^2\mid \Phi_{28341}(28341)$$
