The question is equivalent to asking whether $b(n):=\frac{a(n)}{a(n-1)}$ are squarefree. The sequence $b(n)$ is listed in [OEIS A231831](https://oeis.org/A231831) and satisfies the recurrence $b(n+1) = b(n)^3 + b(n)^2 - 1$ with $b(1)=3$. Its terms are pairwise coprime, implying that each prime divides at most one term.

Quite similarly to [my treatment](https://mathoverflow.net/q/335636) of Sylvester sequence, I have computed all primes $p$ below $10^{10}$ such that $p\mid b(n)$ for some $n$ (there are $16944$ such primes), and verified that $p^2\nmid b(n)$. The primes $p$ and corresponding indices $n$ are now listed in [OEIS A362250](https://oeis.org/A362250) and [OEIS A362251](https://oeis.org/A362251), respectively.

---

The question can also be analyzed heuristically by considering the map $x\mapsto x^3 + x^2 - 1$ modulo prime $p$ as random. Various kinds of statistics for such mappings and their functional graphs are given by [Flajolet and Odlyzko (1989)](https://hal.inria.fr/inria-00075445). My very rough analysis suggests that the "probability" for prime $p$ to divide some $b(n)$ is proportional to $\frac1{p^{1/2}}$, and the "probability" for $p^2$ to divide some $b(n)$ is proportional to $\frac1{p^{3/2}}$. The latter means that one should expect only a finite number of primes $p$ with $p^2\mid b(n)$, if any exists at all.

---

There is a somewhat similar sequence [OEIS A231830](https://oeis.org/A231830) satisfying $c(n+1) = c(n)^3 - c(n)^2 + 1$ with $c(1)=5$. The above analysis applies to this sequence as well. Confirming that neither prime below $10^{10}$ divides its terms when squared is underway.