We consider the numerical stability problem for fractional delay Schrödinger equations involving a Caputo fractional derivative in time, which is developed by Galerkin finite element method (FEM) in space and fractional linear multistep methods (FLMMs) in time. Through rigorous analyses of the characteristic equation yielded by the Laplace transform, we first present an \(\varvec{\alpha }\) -dependent coefficient criterion to ensure the stability of spatially semidiscrete Galerkin FEM and extend the stability property to all convergent spatially semidiscrete methods. Secondly, by introducing a decoupling technique after \(\varvec{\mathcal {Z}}\) transform, we prove the stability of FLMMs generated by both A-stable and A( \(\varvec{\beta }\) )-stable linear multistep methods, without any restriction on the time step size. The stability results are formulated by the fractional exponent, the principal eigenvalue of Dirichlet Laplacian, and the mesh size, but are not related to the delay and time step size. However, for a general spatial region, the principal eigenvalue of Dirichlet Laplacian is always unavailable. In order to provide an effective method for stability detection, when the stability condition is violated, we prove that the fractional trapezoidal rule is an effective method to detect stability because it can not only maintain the stable behavior of the semidiscrete solution, but also the unstable behavior. Extensive numerical experiments for fractional delay Schrödinger equations confirm the long-time decay behaviors of the fully discrete numerical solutions.
