Here's an
observation
I posted on Twitter together with a
proof
.
The Twitter version has some errors which I have fixed, but further corrections are welcome!
Let $M$ be an integer $k\times k$ matrix with distinct eigenvalues $a_{1}, \ldots, a_{k}$, and characteristic polynomial $P(x) = (x-a_{1})\cdots(x-a_{k}) = x^k - \sum_{i=0}^{k-1} b_i x^i$.
Then the equation $f_{n} := \trace\left( M^n \cdot {P'(M)}^{-1}\right) \quad (n = 0,1,2,\ldots)$ defines a sequence of integers.
You can recover the Multinacci numbers (k−ナッチ数) by plugging in a companion matrix for $M$.
For instance if you choose $M = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}$, then
$P(x) = x^{2} -x -1 = (x-\frac{1+\sqrt{5}}{2})(x-\frac{1-\sqrt{5}}{2})$ and $P'(M) = 2M-I = \begin{pmatrix} -1 & 2 \\ 2 & 1 \end{pmatrix}$.
The $f_{n}$ turn out to be the Fibonacci numbers, but ${P'(M)}^{-1} = \frac{1}{5}\begin{pmatrix} -1 & 2 \\ 2 & 1 \end{pmatrix}$ is not an integer matrix.
It is not obvious why the trace should be an integer in general.
The first $k$ values are $f_{0} = \ldots = f_{k-2} = 0$ and $f_{k-1} = 1$.
The proof idea is thanks to @apu_yokai, who linked to his older article about partial fraction decomposition.
As $M$ has $k$ distinct eigenvalues $a_{1},...,a_{k}$ we can act (for the purpose of computing the trace) as if $M$ were diagonal:
$$ f_n = \sum_{i=1}^k \frac{a_i^n}{P'(a_i)} \qquad (n = 0,1,2,\ldots)$$
By partial fraction decomposition we get $f_{0}$:
$$ \sum_{i=1}^k \frac{x}{x-a_i} \frac{1}{P'(a_i)} = \frac{x}{(x-a_{1})\cdots(x-a_{k})} = \frac{x}{P(x)} $$
Take the limit $x \to \infty$:
$$ f_0 = \sum_{i=1}^k \frac{1}{P'(a_i)} = \begin{cases} 0 & \text{ for } 1 < k \\ 1 & \text{ for } 1 = k \\ \end{cases} $$
$f_{1},...,f_{k-1}$ we get inductively.
$n-1 \to n$, while $n < k$:
$$ \underset{=0}{\underbrace{(x^{n-1} f_0 + x^{n-2} f_1 + ... + x f_{n-2} + f_{n-1})x}} + \sum_{i=1}^k \frac{a_i^n x}{x-a_i} \frac{1}{P'(a_i)} $$
$$ = \sum_{i=1}^k \left( (x^{n-1} + x^{n-2}a_i + ... + xa_i^{n-2} + a_i^{n-1})x + \frac{a_i^{n} x}{x-a_i} \right) \frac{1}{P'(a_i)} $$
$$ = \sum_{i=1}^k \frac{(x^{n} - a_i^n)x + a_i^n x}{x-a_i} \frac{1}{\prod_{j\neq i} (a_i-a_j)} = \frac{x^{n+1}}{P(x)} $$
Therefore
$$ f_n = \lim_{x\to\infty} \sum_{i=1}^k \frac{a_i^n x}{x-a_i} \frac{1}{P'(a_i)} = \lim_{x\to\infty} \frac{x^{n+1}}{P(x)} = \begin{cases} 0 & \text{ if } n < k - 1 \\ 1 & \text{ if } n = k - 1 \\ \end{cases} $$
The sequence ${(f_{n})}_{n\geq 0}$ satisfies the recursion scheme $f_{n+k} = \sum_{i=0}^{k-1} b_{i} f_{n+i}$ for all $n\geq 0$.
This is due to Cayley-Hamilton and linearity of the trace
\begin{align} & M^{n}P(M)P'(M)^{-1} = 0 \\ \implies & M^{n+k}P'(M)^{-1} = \sum_{i=0}^{k-1} b_{i} M^{n+i}P'(M)^{-1} \\ \implies & \trace(M^{n+k}P'(M)^{-1}) = \sum_{i=0}^{k-1} b_{i} \trace(M^{n+i}P'(M)^{-1}) \\ \iff & f _{n+k} = \sum _{i=0}^{k-1} b _{i} f _{n+i} \\ \end{align}
Because the initial values are integers and above recursion scheme has integer coefficients all $f_{n}$ are integers.
The condition that $M$ has distinct eigenvalues means that $\gcd(P, P') = 1$, so $P'(M)$ is invertible, and the sequence is well-defined.
Via Lemma 2 the sequence satisfies a recursion scheme with integral coefficients $b_{i}$, and the initial values are also integral by Lemma 1. Therefore all $f_{n}$ are integral by induction.