In this paper we consider higher-order Riesz transforms associated with the Schrödinger operator \(L=-\Delta +P\) on \({\mathbb {R}}^n\) , where the potential P is an arbitrary nonnegative polynomial. We show that for every \(\alpha ,\beta \in {\mathbb {N}}_0^n\) , the transform \(\partial ^\alpha L^{-(|\alpha |+|\beta |)/2}\partial ^{\beta }\) is bounded from \(H_L^p({\mathbb {R}}^n)\) to \(L^p({\mathbb {R}}^n)\) for \(0< p \le 1\) , and if \(|\alpha | \ne 0\) then \(\partial ^\alpha L^{-(|\alpha |+|\beta |)/2}\partial ^{\beta }\) is also bounded from \(H_L^p({\mathbb {R}}^n)\) to \(H^p({\mathbb {R}}^n)\) for \(\frac{n}{n+1}< p \le 1\) , where \(H_L^p({\mathbb {R}}^n)\) and \(H^p({\mathbb {R}}^n)\) are the Hardy space associated with L and the classical Hardy space, respectively.