Generalized parameter-shift rule

Given a Hermitian observable $O$, a Hermitian generator $G$ with eigenvalues $\{\omega_j\}_{j\in[d]}$, the single-parameter gate $U(x) = e^{ixG}$, and a pure state $|\psi\rangle$, consider the loss $\ell(x) = \langle\psi| U^\dagger(x) O U(x) |\psi\rangle$. Let $\{\Omega_p\}_{p\in[R]} := \{\omega_k - \omega_j : j,k\in[d], \omega_k > \omega_j\}$ be the $R$ unique positive eigenvalue differences of $G$. Assuming these frequencies are equidistant integers $\Omega_p = p$ (without loss of generality, by rescaling), the first and second derivatives of $\ell$ at the origin are given by the generalized parameter-shift rule $$\ell'(0) = \sum_{\mu=1}^{2R} \ell\left(\frac{2\mu-1}{2R}\pi\right) \frac{(-1)^{\mu-1}}{4R \sin^2\left(\frac{2\mu-1}{4R}\pi\right)},$$ $$\ell''(0) = -\ell(0)\frac{2R^2+1}{6} + \sum_{\mu=1}^{2R-1} \ell\left(\frac{\mu\pi}{R}\right) \frac{(-1)^{\mu-1}}{2\sin^2\left(\frac{\mu\pi}{2R}\right)}.$$ For $R=1$ and $R=2$ these reduce to the standard two-term and four-term parameter-shift rules.