Let $A\in \mathsf{GL}(4,\mathbb{R})$ be the following matrix:
$$A=\begin{pmatrix} \cos a&-\sin a&0&0\\ \sin a&\cos a&0&0\\0&0&\cos b&-\sin b\\0&0&\sin b&\cos b \end{pmatrix}$$
Assume that $A$ has finite order, i.e. $a$ and $b$ are rational multiples of $\pi$.
The problem is to give necessary conditions on the characteristic polynomial of $A$ to make it integer, i.e. $P_A(\lambda)\in\mathbb{Z}[\lambda]$
The characteristic polynomial is $P_A(\lambda)=\lambda^4-2\lambda^3(\cos a+\cos b)+\lambda^2(2+4\cos a\cos b)-2\lambda(\cos a+\cos b)+1$ (it is a symmetric polynomial).
Now we want that
\begin{cases} 2(\cos a+\cos b)\in\mathbb{Z}\\ 4\cos a\cos b\in \mathbb{Z} \end{cases}
Question: In this case I've solved the system by hand and found the possible values of $a$ and $b$ but I wonder if there is a more efficient way to do this, or if there is some deeper theory involved, because I want to generalize to higher values (for example when there are three blocks, $a$, $b$, $c$ and so on). I see that in someway elementary symmetric polynomials appear but I don't know if this helps.
Any comment or suggestion will be appreciated!
We'll say that $A$ is block diagonal with $k$ blocks of size $2 \times 2$ with each $2 \times 2$ block in the sine/cosine form you present, and $A$ has finite order.
The problem of counting the possible characteristic polynomials of $A$ amounts to counting the polynomials $p(x)$ of degree $2k$ for which $p(x) \mid x^n - 1$ for some positive integer $n$ and $p(x)$ has no real roots. Equivalently, we want the number of degree $2k$ polynomials that can be written as a product of cyclotomic polynomials with degree at least $2$.
Because the degree the cyclotomic polynomial $\Phi_n$ is $\varphi(n)$, where $\varphi$ denotes the totient function, we can reframe the problem as follows: we want to count the number of multisets (sets with repetition) $S \subset \{3,4,5,\dots\}$ that satisfy the condition $$ \sum_{j \in S} \varphi(j) = 2k. $$ Because the totient function has lower bounds such as $\varphi(n) \geq \sqrt{n/2}$, this problem could be solved by "brute force". For example, it is sufficient to take $S \subset \{3,4,\dots, 8k^2\}$, and the number of elements in $S$ is at most $k$.
This is also made much easier if we happen to know the number of integers satisfying $\varphi(n) = j$ for $1 \leq j \leq 2k$, as in the list given here or in this OEIS entry.
To address your $k = 2$ problem in detail, it suffices to note that there are $3$ choices of $n$ with $\varphi(n) = 2$ (namely $3,4,6$), and $4$ choices of $n$ with $\varphi(n) = 4$ (namely $5,8,10,12$). Thus, the total number of possible characteristic polynomials for the matrix satisfying the criteria of your problem will simply be $$ \binom{3 + 2 - 1}{2} + 4 = 6 + 4 = 10, $$ where we make use of the multiset coefficient here. The corresponding number of matrices is $$ 3^2 + 4 = 13. $$