I am an undergrad math major and have been doing some research as a hobby. In short, I am looking at primes whose squares divide some number of the form $n^n + (-1)^n (n-1)^{(n-1)}$. I have found that for primes up to 800,000, about 55% of them are equivalent to 1 mod 4. I know there are theorems that prove that primes in general are evenly distributed across any modulus; more specifically, that the density of primes equivalent to a mod n is $\dfrac{1}{\phi(n)}$. I thought it was very odd that such a small but seemingly solid bias would appear in so random a place - does anyone know of similar subsets of primes where it has been proven/conjectured that some form of bias like this appears?
As requested by a comment, I am adding some information about these primes here (A283146 in the OEIS). As a means of ease of notation, let $X(n) = n^n + (-1)^n(n-1)^{n-1}$.
For a given prime $p$, if $p^2 | X(n)$, then $p^2 | X(a)$ for all $a \equiv n \pmod{p(p-1)}$. A similar statement holds for non-divisibility, so that all solutions can be determiend by testing each residue class mod $p(p-1)$. Without listing the entire proof, this can be proven by simplifying $X(n+p(p-1)$ using the binomial theorem and using the Euler-Fermat Theorem to eliminate some coefficients. These simplifications show that $X(n+p(p-1)) \equiv X(n) \pmod{p(p-1)}$ independent of whether $p^2 | X(n)$.
If $p^2 | X(n)$, then $p^2 | X(p(p-1) + 1 - n)$. That is, residue classes $a,b$ mod $p(p-1)$ that sum to 1 mod $p(p-1)$ either $p^2$ divides both $X(a)$ and $X(b)$ or neither. The proof for this statement uses a very similar approach to the first proof; it uses the binomial theorem to simplify $X(p(p-1) + 1 - n)$ and uses Euler-Fermat to deal with some leftover coefficients.
There is a correspondence between consecutive $p$th powers mod $p^2$ and the values of $n$ for which $p^2 | X(n)$. It should be noted that the polynomial $f_p(x) = \dfrac{(x+1)^p - x^p - 1}{p} \mod{p}$ detects these roots. There is a generalization of this statement in the paper referenced in the OEIS article. This paper proves this statement by explicit construction of bijections between these, though I honestly am not enough of an expert to detail the proof here.
It has been conjectured that the relative density of this set of primes within the primes is around 15%. This has followed from numerical evidence and a heuristic about the roots of $f_p(x)$ that I do not fully understand but are detailed in the paper I cite.
If we let $\mathcal{C}_p = \{ n \in \mathbb{Z}_{p(p-1)}$ so that $p^2 | X(n) \}$, it is conjectured that $|\mathcal{C}_p|$ averages to 1 over the primes. This is a heuristic used in the cited paper and is explained using a probabalistic argument.
Any help would be much appreciated! Here is the link to the paper I reference