We study the incompressible magnetohydrodynamic system endowed with a Caputo time-fractional derivative of order \(0<\alpha <1\) . First, we reformulate the coupled momentum-induction equations by applying a double-curl projection that eliminates both the hydrodynamic pressure and the magnetic pseudo-pressure, producing a divergence-free velocity-magnetic field pair. Next, we then discretise the reformulated problem by combining a divergence-free Fourier spectral approximation in space with a variable-step L1 convolution discretisation in time. Our fully discrete method satisfies a discrete fractional kinetic–magnetic energy inequality, enforces the divergence constraints to machine precision, and reduces to the classical magnetohydrodynamic energy law as \(\alpha \rightarrow 1^{-}\) , thereby ensuring asymptotic compatibility. Using discrete orthogonal-complementary convolution identities together with discrete Sobolev embeddings, we derive optimal error bounds of order \(O\!\bigl (\tau ^{2-\alpha }+N^{-m}\bigr )\) on arbitrarily graded time meshes. Finally, we present numerical experiments, including fractional magnetic Taylor–Green and Orszag–Tang vortices, which confirm the theoretical convergence rates, demonstrate monotone energy decay, and highlight the efficiency of the adaptive time-stepping strategy.
