In this paper, we present unconditionally optimal error estimates of linearized Crank–Nicolson Galerkin finite element methods for a strongly nonlinear parabolic system in $$\mathbb {R}^d\ (d=2,3)$$ . However, all previous works required certain time-step conditions that were dependent on the spatial mesh size. In order to overcome several entitative difficulties caused by the strong nonlinearity of the system, the proof takes two steps. First, by using a temporal-spatial error splitting argument and a new technique, optimal $$L^2$$ error estimates of the numerical schemes can be obtained under the condition $$\tau \ge h$$ , where $$\tau $$ denotes the time-step size and h is the spatial mesh size. Second, we obtain the boundedness of numerical solutions by mathematical induction and inverse inequality when $$\tau \le h$$ . Then, optimal $$L^2$$ and $$H^1$$ error estimates are proved in a different way for such case. Numerical results are given to illustrate our theoretical analyses.