Let $A$ be an integer valued matrix, with $\det{A} \neq 0$. Show that there exists an infinite number of positive integers $p$ for which there exists a constant $C_p$ such that all entries of the matrix $A^{C_p} - I_2$ are divisible by $p$, where $I_2$ is the second order identity matrix.
This is what I have been able to prove:
1) If the matrix is diagonal, with diagonal elements, say $a$ and $b$, then by choosing $p$ coprime to $a$ and $b$ and choosing $C_p = \phi(p)$, all elements of $A^{C_p} - I_2$ are divisible by $p$, due to Euler's theorem.
2) If the matrix is, say, $$ A = \begin{bmatrix} a & b\\ 0&c \end{bmatrix},$$ then if $a \neq c$ (the $a = c$ case is simple to prove) we can show that $$A^n = \begin{bmatrix} a^n & b\left(\frac{a^n - c^n}{a-c}\right) \\ 0 & c^n\end{bmatrix}.$$
Choosing $p$ coprime to $a,c, a-c$, and $C_p = \phi(p)$, then the problem is again solved by Euler's theorem.
I cannot prove it if all entries are nonzero. I have been able to show that the gcd of the entries of $A - I_2$ is $\leq $ the gcd of the entries of $A^2 - I_2$. If I could prove that the sequence: $$a_n = \text{gcd of entries of } A^{2^n} - I_2$$ does not eventually become stationary, then the problem is solved. Experimenting with Mathematica shows that this is the case.
Any ideas?
I will look for some reference for this: given a nonzero integer $n,$ there are infinitely many primes for which $n$ is a square (a quadratic residue). The two proofs i can think of use quadratic reciprocity, Dirichlet's result on primes in arithmetic progressions, or Chebotarev density.
https://number.subwiki.org/wiki/Every_integer_is_a_quadratic_residue_for_infinitely_many_primes
they use:
https://number.subwiki.org/wiki/Set_of_prime_divisors_of_values_of_nonconstant_polynomial_with_integer_coefficients_is_infinite
Given a matrix $$ \left( \begin{array}{cc} \alpha & \beta \\ \gamma & \delta \end{array} \right) $$ the integer in question is $$ n = (\alpha - \delta)^2 + 4 \beta \gamma $$
If we take an odd prime $p$ that does not divide the determinant $\alpha \delta - \beta \gamma$ or $n,$ and for which $n$ is a quadratic residue, then the eigenvalues of the matrix lie in the field with $p$ elements as $$ \frac{\alpha + \delta \pm \sqrt n}{2}. $$ Then we can diagonalize the matrix as some $P^{-1}DP = A$ over the field with $p$ elements. There is an exponent such that $D^C = I \pmod p.$ The same applies to $A.$