  The paper focuses on unconditionally optimal error analysis of the fully
discrete Galerkin finite element methods for a general nonlinear parabolic
system in $\R^d$ with $d=2,3$. In terms of a corresponding time-discrete system
of PDEs as proposed in \cite{LS1}, we split the error function into two parts,
one from the temporal discretization and one the spatial discretization. We
prove that the latter is $\tau$-independent and the numerical solution is
bounded in the $L^{\infty}$ and $W^{1,\infty}$ norms by the inverse
inequalities. With the boundedness of the numerical solution, optimal error
estimates can be obtained unconditionally in a routine way. Several numerical
examples in two and three dimensional spaces are given to support our
theoretical analysis.
