Let m be a nonnegative integer, and let \(n\ge 2^{m+1}+1.\) In this paper, we consider the higher order Schrödinger type operator \({\mathcal {H}}_{2^m}=(-\Delta )^{2^m}+V^{2^m} \) on \({\mathbb {R}}^n,\) and establish the \(L^p({\mathbb {R}}^n)\) boundedness of Riesz transforms \(\nabla ^j {\mathcal {H}}_{2^m}^{-\frac{j}{2^{m+1}}} (j=1,2,\cdot \cdot \cdot ,2^{m+1}-1)\) and their commutators. Here, V is a nonnegative potential belonging to both the reverse Hölder class \(RH_s\) for \(s \ge \frac{n}{2}\) , and the Gaussian class associated with \((-\Delta )^{2^m}\) .