Prime Appearances in Fibonacci Number Factorizations
A partial result: any FPA mus be congruent to $2, 3 \bmod 5$, and these primes have density $\frac{1}{2}$ in the set of all primes by the strong form of Dirichlet's theorem.
First, I claim that Binet's formula $F_n = \frac{\phi^n - \varphi^n}{\phi - \varphi}$ (my numbering differs from yours), where $\phi, \varphi$ are the two roots of $x^2 = x + 1$, is valid modulo $p$ for any prime $p \neq 5$, where $\phi, \varphi$ are interpreted as elements of the splitting field of $x^2 = x + 1$ over $\mathbb{F}_p$. The proof is the same proof as in the usual case, the point being that the discriminant of $x^2 = x + 1$ is $5$, so for any other prime there are two distinct roots. Moreover, since $(\phi - \varphi)^2 = 5$, we know that $\phi - \varphi \neq 0$ for any such prime. It follows that $F_n \equiv 0 \bmod p$ if and only if $\phi^n \equiv \varphi^n \bmod p$.
If $\left( \frac{5}{p} \right) = 1$, then $x^2 = x + 1$ splits over $\mathbb{F}_p$, from which it follows that $\phi^{p-1} \equiv \varphi^{p-1} \equiv 1 \bmod p$, hence that $$F_{p-1} \equiv 0 \bmod p.$$
By quadratic reciprocity, $\left( \frac{5}{p} \right) = \left( \frac{p}{5} \right)$, hence these are precisely the primes $p \equiv 1, 4 \bmod 5$. So none of these primes are FPAs.
If $\left( \frac{5}{p} \right) = -1$, then $x^2 = x + 1$ has splitting field $\mathbb{F}_{p^2}$. The Galois group is generated by the Frobenius map $x \mapsto x^p$, hence $\phi^p \equiv \varphi \bmod p$ and, similarly, $\varphi^p \equiv \phi \bmod p$, hence $\phi^{p+1} \equiv \phi \varphi \equiv \varphi^{p+1} \equiv -1$, from which it follows that $$F_{p+1} \equiv 0 \bmod p.$$
This occurs if and only if $p \equiv 2, 3 \bmod 5$.
Edit: Second partial result: any FPA must be congruent to $3 \bmod 4$. Now the density has been reduced to $\frac{1}{4}$. To see this, note that since $\phi \varphi = -1$, the condition that $\phi^n \equiv \varphi^n \bmod p$ is equivalent to the condition that $(-\phi^2)^n \equiv 1 \bmod p$. On the other hand, we know that when $p \equiv 2, 3 \bmod 5$ we have $\phi^{p+1} \equiv -1 \bmod p$, hence $$\left( - \phi^2 \right)^{ \frac{p+1}{2} } \equiv (-1)^{ \frac{p-1}{2} } \bmod p.$$
It follows that when $p \equiv 1 \bmod 4$ we have
$$F_{ \frac{p+1}{2} } \equiv 0 \bmod p.$$
I don't expect these necessary conditions to be sufficient. The question is similar to asking that $- \phi^2$ behave like a primitive root, so this question or a variation on it may end up being an open problem.
Edit #2: For example, you can verify for yourself that $F_{16} \equiv 0 \bmod 47$ even though $47 \equiv 2 \bmod 5$ and $47 \equiv 3 \bmod 4$.
Building on Qiaochu Yuan's answer, the problem is not to check that $\alpha = -\phi^2$ behaves like a primitive root, but to check that it has no more $n$th roots in $\mathbb{F}_{p^2}$, other than the ones given by its $(p-1)$th root.
Actually, you have already proved the partial result $\alpha$ has exactly enough nested square roots if and only if $p = 3 \mod 4$, so with probability $\frac{1}{2}$.
If $n$ is not a multiple of $2$ or $5$, pick a primitive $n$th root of unity $\zeta_n$, and we can look at the extension $\mathbb{Q} \subset \mathbb{Q}(\alpha,\zeta_n,\beta = \sqrt[n]\alpha)$
For any prime not dividing $n$, and different from $2$ or $5$, denote by $\sigma$ the Frobenius morphism. Chebotarev's theorem implies that the density of the primes for any possible action of $\sigma$ is always the same ($\frac1{2n\varphi(n)}$) :
$\sigma(\alpha) = \alpha^\epsilon$ with $\epsilon = \pm 1$
$\sigma(\zeta_n) = \zeta_n^{p \mod n}$ (this is Dirichlet's theorem)
$\sigma(\beta) = \zeta_n^k \beta^\epsilon : \sigma(\sigma(\beta)\beta^{-\epsilon})^n = \sigma(\alpha^{\epsilon-\epsilon}) = 1$, so $\sigma(\beta)\beta^{-\epsilon} = \zeta_n^k $ for some $k \in \mathbb{Z}/n\mathbb{Z}$.
Then for any divisor $d$ of $n$, the $d$th roots of $\alpha$ in $\mathbb{F}_{p^2}$ are the $\rho_i = (\zeta^i\beta)^{n/d}$ such that $\rho_i = \sigma\sigma(\rho_i) $. But $$ \sigma\sigma(\rho_i^{n/d}) = \sigma\sigma(\zeta_n^{in/d}\beta^{n/d}) =\sigma(\zeta_n^{(pi+k)n/d}\beta^{-n/d}) = \zeta_n^{(p^2i+(p-1)k)n/d}\beta^{n/d} = \zeta_n^{(p-1)((p+1)i+k)n/d} \rho_i$$ So it was equivalent to $(p-1)((p+1)i+k)n/d = 0 \mod n$, which is also $(p-1)((p+1)i+k) = 0 \mod d$ : The number of $d$th roots in the number of classes $i$ mod $d$ such that $(p-1)((p+1)i+k) = 0 \mod d$
This agrees with what we know, for example when $d$ is prime, if $p = 1 \mod d$, then since $\alpha$ has $(p-1)$th roots, it has $d$th roots ; if $p \neq \pm 1 \mod d$, then $x \mapsto x^d$ is a permutation of $\mathbb{F}_{p^2}$ so we can always find exactly one root ; and if $p = -1 \mod d$ then it has $d$th roots if and only if $k=0 \mod d$, which happens 1/d of the time.
The point of having a generic $n$ instead of a prime is that it tells us having $q$th roots is "independant" of having a $q'$th roots for coprime $q$ and $q'$.
Now we can say that for primes $p$ other than $2$ or $5$:
- with probability $\frac12$, $\alpha \notin \mathbb{F}_{p^2}$ and $p$ is therefore disqualified from being a PFA.
The other $\frac{1}{2}$ is split independantly for all primes $q$ other than $2$ or $5$ :
- with probability $\frac1{(q-1)}$, if $p = 1 \mod q$ and $\alpha$ has $q$th roots in $\mathbb{F}_{p^2}$.
- with probability $1-\frac2{q-1}$, $p^2 \neq \pm 1 \mod q$ and $\alpha$ has one $q$th root in $\mathbb{F}_{p^2}$.
- with probability $\frac1q$, $p = -1 \mod q$ and $\alpha$ has no $q$th root in $\mathbb{F}_{p^2}$.
with probability $\frac1{q(q-1)}$, $p = -1 \mod q$ and $\alpha$ has $q$th roots in $\mathbb{F}_{p^2}$, and $p$ is therefore disqualified from being a PFA.
for $q=5$, there is no obstruction to being a PFA
- for $q=2$, $p$ is not a PFA with probability $\frac12$
Thus I expect the density of PFA to be precisely $\frac{10}{19} \Pi (1-\frac1{p(p-1)})$ where $p$ ranges over all primes.