In this paper, we examine and provide numerical solutions to the nonlinear fractional order time-space diffusion equations with the influence of temporal delay. An effective high-order numerical scheme that mixes the so-called Alikhanov L2 − 1σ formula side by side to the power of the Galerkin method is presented. Specifically, the time-fractional component is estimated using the uniform L2−1σ difference formula, while the spatial fractional operator is approximated using the Legendre-Galerkin spectral approximation. In addition, Taylor’s approximations are used to discretize the term of the nonlinear source function. It has been shown theoretically that the suggested scheme’s numerical solution is unconditionally stable, with a second-order time-convergence and a space-convergent order of exponential rate. Furthermore, a suitable discrete fractional Grönwall inequality is then utilized to quantify error estimates for the derived solution. Finally, we provide a numerical test that closely matches the theoretical investigation to assess the efficacy of the suggested method.