I have a problem wherein I have defined a function $I_r(t) = \int e^{(2r-1)at} \int e^{(2r-3)at} \cdots \int e^{at} dt\cdots dt$, and $I_r(0) = 0$, for $r = 1,2,3,\ldots$.

I find that $e^{-ar^2t} I_r(t) = \left(1-e^{-at}\right)^r q(t)$, where $q(exp(-at))$ is a polynomial of $e^{-at}$. Is there a general technique for evaluating repeated integrals of this type that allows me to write $q$ in a nice clean way?

If I took $I^*_r(t) = \int e^{at} \int e^{at}\cdots \int e^{at} dt\cdots dt$ with $I^*_r(0) = 0$, and multiplied by $e^{-art}$, I would get $e^{-art}I^*_r(t) = (1-e^{-at})^r$. I am looking for a nice closed-form solution where I have a quadratic in $r$.

This is related to the derivation of a discrete probability distribution where the transition rate function is quadratic w/r.t. the number of events per cell.