For any $N\times N$ matrix $A$ we set $\DeclareMathOperator{\tr}{tr}$ $\newcommand{\bE}{\mathbb{E}}$

$$ |A|^2:=\sum_{i,j}|a_{ij}|^2 =\tr(A A^*). $$

There exists a constant $C=C(N)>0$ such that for any complex $N\times N$ matrices $A, B$ we have

$$ |AB|\leq C |A|\cdot |B|. $$

In particular, if $A_1,\dotsc, A_n$ are complex $N\times N$ matrices we have

$$ |A_1\cdots A_n|\leq C^{n-1}|A_1|\cdots |A_n|. $$

We have

$$ \Bigl\vert \bE(AX)^n\Bigr\vert\leq \bE\bigl(\;|(A(X)^n|\;\bigr) \leq C^{2n-1}|A|^n \bE\bigl(\;|X|^n\;\bigr) =C^{2n-1}|A|^n \bE\bigl(\;|Z|^{2n}\;\bigr) $$

$$ = C^{2n-1}|A|^n\bE\Bigl( \;\Bigl(\; \sum_{i,j}|z_{ij}|^2\;\Bigr)^n\;\Bigr). $$

Now the question reduces to the following. Suppose we are given independent, complex Gaussian variables $\zeta_1,\dotsc, \zeta_m$. Are the moments of $|\zeta_1|^2+\cdots +|\zeta_m|^2$ finite? The answer is yes, because the probablity

$$ p(R)=P(|\zeta_1|^2+\cdots +|\zeta_m|^2 > R^2) $$

goes to zero exponentially fast since the random vector $(\zeta_1,\dotsc,\zeta_m)$ is Gaussian.

Next use the fact that for any nonnegative random variable $X$ we have the equality

$$ \bE(X^k) =\int_0^\infty kx^{k-1} P(X>x) dx,\;\;\forall k\geq 1. $$