What is the structure of invariant matrix polynomials?

Yes. This is almost true in general, but not quite for two reasons: one involving the size of the field and the other involving characteristic.

The correct statement over an infinite field $F$ is that the algebra of invariant polynomials is generated by the elementary symmetric functions of the eigenvalues, or more concretely by the coefficients of the characteristic polynomial $\det (I - At)$. To see this, pass WLOG to the algebraic closure $\bar{F}$ (which we can do because $P$ is polynomial). The identity $P(gAg^{-1}) = P(A)$ is a polynomial identity in the coefficients of $A$, the coefficients of $g$, and $\det g$, and since $F$ is an infinite field two polynomials which are equal when all possible elements of $F$ are plugged in are equal identically. Consequently, this identity continues to hold over $\bar{F}$. (This part of the argument can be ignored over $\mathbb{C}$.)

Now suppose that $M$ is a diagonalizable matrix. Then $P(M)$, being conjugation-invariant, is necessarily a polynomial, in fact a symmetric polynomial, in the eigenvalues of $M$. Since the diagonalizable matrices over $\bar{F}$ are Zariski-dense (as they contain the matrices with distinct eigenvalues, which is a Zariski-open condition), it follows that $P$ must be a symmetric polynomial in the eigenvalues of $M$ identically. Hence $P$ is a polynomial in the elementary symmetric polynomials with coefficients in $\bar{F}$, but using the fact that $P$ is defined over $F$ and letting $M$ be a suitable collection of companion matrices shows that $P$ is a polynomial in the elementary symmetric polynomials with coefficients in $F$.

The above argument does not work over a finite field because the identity $P(gAg^{-1}) = P(A)$ is no longer guaranteed to continue to hold over the algebraic closure, although I don't know a counterexample in this setting.

The result you want is true in characteristic zero or characteristic greater than $n$ by Newton's identities but false in small positive characteristic because using Newton's identities will require dividing by zero. For example, in characteristic $2$, unlike in higher characteristic, it is not possible to express $$e_2 = \sum_{i < j} \lambda_i \lambda_j$$

in terms of $$p_1 = \sum_i \lambda_i$$

and $$p_2 = \sum_i \lambda_i^2 = (\sum_i \lambda_i)^2 = p_1^2$$

(where $\lambda_i$ are the eigenvalues).