I would like a closed form solution to $$ P(n):=\sum_{m=0}^\infty \textrm{NegBin}_{r,p}(m)\textrm{Poisson}_{m\lambda}(n)\\ = \sum_{m=0}^\infty \left(\begin{array}{c}{m+r-1} \\ {m}\end{array}\right) \cdot(1-p)^{r} p^{m}\frac{(\lambda m)^{n} e^{-\lambda m}}{n !} $$ i.e. a Poisson distribution for random variable $n$ whose scale parameter is proportional to a negative binomially distributed random variable, $m$, that is marginalized over. $\lambda>0$ is the proportionality constant. $r$ and $p$ are the standard parameters of the negative binomial distribution.
Given that a convolution of negative binomials has a closed form (https://www.tandfonline.com/doi/full/10.1080/03610926.2015.1053931) I feel like the sum should be doable, I just don't know through what method.
The gamma function representation of the binomial coefficient in the negative binomial, $\frac{\Gamma(r+m)}{m ! \Gamma(r)}$, might be a way, but that seems like overkill...maybe is there just some summation trick, but so far I haven't found it.
The hierarchical model is: $$N \mid M \sim \operatorname{Poisson}(\lambda M), \quad M \sim \operatorname{NegBinomial}(r,p)$$ where $\lambda, r, p$ are fixed parameters, and $M$ is parametrized such that $M \in \{0, 1, 2, \ldots \}$. Then the PMF of the unconditional distribution of $N$ is $$\Pr[N = n] = \sum_{m=0}^\infty \Pr[N = n \mid M = m]\Pr[M = m] = \sum_{m=0}^\infty e^{-\lambda m} \frac{(\lambda m)^n}{n!} \binom{m+r-1}{r-1} p^r (1-p)^m .$$ (Note that your binomial coefficient is equivalent to mine.) Then pulling out the constant terms with respect to the index of summation, $$\Pr[N = n] = \frac{\lambda^n p^r}{n!} \sum_{m=0}^\infty m^n \binom{m+r-1}{r-1} (e^{-\lambda} (1-p))^m.$$ This suggests letting $$p^* = 1 - e^{-\lambda}(1-p),$$ so that $$\Pr[N = n] = \frac{\lambda^n}{n!} \left(\frac{p}{p^*}\right)^r \sum_{m=0}^\infty m^n \binom{m+r-1}{r-1} (p^*)^r (1 - p^*)^m = \frac{\lambda^n}{n!} \left(\frac{p}{p^*}\right)^r \operatorname{E}[(M^*)^n],$$ where $$M^* \sim \operatorname{NegBinomial}(r,p^*).$$ From the MGF of the negative binomial distribution, we can write this expectation as $$\operatorname{E}[(M^*)^n] = \mathbb M_{M^*}^{(n)}(0) = \frac{d^n}{dt^n} \left[ \left(\frac{p^*}{1-e^t(1-p^*)}\right)^r \right]_{t=0}.$$ These derivatives are quite complicated, unfortunately.