Abstract
We study the incompressible magnetohydrodynamic system endowed with a Caputo time-fractional derivative of order [Formula: see text]. 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 [Formula: see text], thereby ensuring asymptotic compatibility. Using discrete orthogonal-complementary convolution identities together with discrete Sobolev embeddings, we derive optimal error bounds of order [Formula: see text] 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.