Integrals involving derivatives of Legendre polynomials frequently arise in applications ranging from multipole expansions for processes involving electromagnetic probes to spectral methods in numerical physics. Despite their practical relevance, closed-form expressions for such integrals - particularly involving arbitrary derivative orders - are not readily accessible in standard references or symbolic tools. In this note, we derive and present general analytic expressions for integrals of the form -₁^+1 dx P^ (q) ₍ (x) P^ (k) ₌ (x), where P₍ (x) and P₌ (x) are Legendre polynomials and q, k denote their order of differentiation. Using repeated integration by parts, parity arguments, and closed-form boundary evaluations, we obtain explicit binomial and Gamma-function representations valid for all non-negative integers n, m, q, k. These results unify and extend known orthogonality relations and provide ready-to-use tools for analytic and computational contexts.
Wunderlich et al. (Fri,) studied this question.