Abstract: Direction of Arrival (DOA) estimation is a core technique in array signal processing. Although sparse recovery-based DOA estimation methods using programmable metasurfaces can effectively reduce system cost and complexity, their performance is affected by the grid mismatch caused by spatial discretization. To address this issue, this paper proposes the off-grid DOA estimation methods for programmable metasurfaces based on joint sparse recovery. Firstly, a joint overcomplete dictionary containing both grid-term and derivative-term is constructed based on the Taylor expansion, thereby establishing a one-dimensional joint DOA estimation model. Building upon this, the model is extended to a two-dimensional scenario, achieving off-grid estimation of both azimuth and elevation angles. To reduce the computational complexity of the two-dimensional model, a model simplification method based on Kronecker decomposition is further introduced, transforming high-dimensional operations into lower-dimensional ones. To solve the formulated models, this paper proposes two joint sparse recovery algorithms, namely Joint Orthogonal Matching Pursuit and Joint L₁-norm minimization. Simulation and practical measurement results demonstrate that the proposed method improves the DOA estimation accuracy of metasurfaces compared to existing typical methods. Meanwhile, the simplified model based on Kronecker decomposition reduces computational complexity while maintaining accuracy comparable to the original model.