In this paper, we study the Orlicz estimates for the parabolic Schrödinger operator

$$L={\mathrm{\partial }}_t-{\mathrm{\Delta }}_X+V,$$

where the nonnegative potential $V$ belongs to a reverse Hölder class on nilpotent Lie groups $\mathbb{G}$ and ${\mathrm{\Delta }}_X$ is the sub-Laplace operator on $\mathbb{G}$. Under appropriate growth conditions of the Young function, we obtain the regularity estimates of the operator $L$ in the Orlicz space by using the domain decomposition method. Our results generalize some existing ones of the $L^p$ estimates.