How were the factors of $\frac{521^{521}-1}{520}$ found?

I believe there is a way to generate the reciprocal polynomial factors $\Phi_n(nx^2)$ or $\Phi_{2n}(nx^2)$ for prime $n$ depending on the congruence $\pmod 4$. The former are $1$ while the latter are $3$ $\pmod 4$. Given the field of $n$-th roots of unity, let $r$ be any primitive $n$-th root of unity, and $f$ is the first factor of $\Phi_n(nx^2)$ or $\Phi_{2n}(nx^2)$. Then there is a polynomial $t$ in terms of $r$ such that the minimal polynomial of $t$ is $f$. It is possible (but may require extensive computation), to generate this minimal polynomial for large prime $n$.

For example, take $n=5$, and

$\Phi_5(5x^2)=(25x^4 - 25x^3 + 15x^2 - 5x + 1)(25x^4 + 25x^3 + 15x^2 + 5x + 1)$

The reciprocal of the first factor is $x^4 - 5x^3 + 15x^2 - 25x + 25$.

If $r$ is a primitive $5$th root of unity, then the minimal polynomial of $r^3-r^2-r+1$ is $x^4 - 5x^3 + 15x^2 - 25x + 25$. When $r$ is a $10$th root of unity, the minimal polynomial $r^3-r^2-r+1$ is $x^4 - 5x^3 + 5x^2 + 5x + 5$. (I thought this would generate the second factor but I guess not). The computations were performed with PARI/GP:

(14:20) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(5)))
%76 = x^4 - 5*x^3 + 15*x^2 - 25*x + 25
(14:20) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(10)))
%77 = x^4 - 5*x^3 + 5*x^2 + 5*x + 5
(14:20) gp >

I also treated $r$ as an $n$-th root of unity for other prime $n$, taking the minimal polynomial of $r^3-r^2-r+1$ and obtained similar polynomials. I don't know weather these polynomials have any relation with the aurifeuillan factors of $\Phi_n(nx^2)$ or $\Phi_{2n}(nx^2)$ or could be used to solve factorization of $(n^n+-1)/(n+-1)$:

(14:22) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(5)))
%78 = x^4 - 5*x^3 + 15*x^2 - 25*x + 25
(14:22) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(7)))
%79 = x^6 - 7*x^5 + 21*x^4 - 35*x^3 + 49*x^2 - 49*x + 49
(14:22) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(11)))
%80 = x^10 - 11*x^9 + 55*x^8 - 165*x^7 + 341*x^6 - 506*x^5 + 484*x^4 - 242*x^3 + 121*x^2 + 121
(14:22) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(13)))
%81 = x^12 - 13*x^11 + 78*x^10 - 286*x^9 + 715*x^8 - 1300*x^7 + 1833*x^6 - 2028*x^5 + 1521*x^4 - 507*x^3 + 169*x^2 + 169*x + 169
(14:22) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(17)))
%82 = x^16 - 17*x^15 + 136*x^14 - 680*x^13 + 2380*x^12 - 6188*x^11 + 12393*x^10 - 19516*x^9 + 24310*x^8 - 24276*x^7 + 20808*x^6 - 16473*x^5 + 10404*x^4 - 2312*x^3 + 1156*x^2 + 1156*x + 289
(14:22) gp > minpoly(Mod(x^3-x^2-x+1,polcyclo(19)))
%83 = x^18 - 19*x^17 + 171*x^16 - 969*x^15 + 3876*x^14 - 11628*x^13 + 27132*x^12 - 50426*x^11 + 75905*x^10 - 93176*x^9 + 92416*x^8 - 73644*x^7 + 50901*x^6 - 35739*x^5 + 23104*x^4 - 4332*x^3 + 3249*x^2 + 2166*x + 361

Any leads further from this?

That formula, $$\frac{n^n - 1}{n - 1},$$ looks so familiar. Maybe it has something to do with my question (with a soon-to-expire bounty) about reciprocating pseudoprimes.

You asked for "any ideas." The first idea I often tell people is to look it up in the OEIS. But first I ran the query (n^n - 1)/(n - 1) on factordb, with slight trepidation that the server could blow up at having to divide by 0 (for $n = 1$).

From that factordb query, I got enough numbers to do an OEIS search and get just one result: http://oeis.org/A023037

For $n \geq 1$, $a(n)$ is the number whose base $n$ representation is a string of $n$ ones. For example, 11111 in base 5 is $a(5) = 781$. - Melvin Peralta, May 23 2016

Of course! $a(n)$ is a base $n$ repunit. And repunits are almost always composite. If $b$ is the base of numeration and $n = b$, then $$\frac{b^n - 1}{n - 1} = \sum_{i = 0}^{n - 1} b^i,$$ and... oops, sorry, that was not as fruitful as I had hoped.

Okay, so we know that this base $n$ repunit is divisible by any factors of $n^n - 1$ that are coprime to $n - 1$. Then, for 521, we see that $520 = 2^3 \times 5 \times 13$.

And then we can confirm that $521^{521} \equiv 1 \pmod 5$, $521^{521} \equiv 1 \pmod 8$ and $521^{521} \equiv 1 \pmod{13}$... sorry, another dead end.

So when is this base $n$ repunit prime? Clearly $n$ itself must be prime. I increased my factordb query to $n = 200$ and got 2, 3, 19, 31, which I feared would give me too many results in the OEIS. It just got me four, and the first one was http://oeis.org/A088790

The next prime is 7547. Still, this does nothing to diminish my astonishment at the discovery of the two factors for the base 521 repunit. There must be something algebraic going on. Maybe something of what I've said helps you, or maybe I've just given you a bunch of red herrings.