Here are my notes on the parameterisation of the Stan negative binomial models, including some discussion of how negative binomial distributions are combined:
negative-binomial.pdf (from here https://github.com/dragonfly-science/negative-binomial). The motivation of this was around figuring out (and checking) how the different parameterisations related to one another.
The short story, is that in the size and probability parameterisation, the sum of draws from two distributions \mathrm{NB}(r, p) and \mathrm{NB}(s, p) is a draw from a distribution \mathrm{NB}(r+s, p). There is a derivation of this here: https://math.stackexchange.com/questions/1054048/negative-binomial-distribution-sum-of-two-random-variables
The parameterisation used in Stan has \mu = r (1-p)/p and \phi = r, so the sum of distributions \mathrm{NB}(\mu_s, \phi_s) and \mathrm{NB}(\mu_r, \phi_r) has \mu = (r+s)(1-p)/p, and \phi=(r+s). If you were adding together N draws with the same \mu and \phi, then the summed distribution would have parameters N\mu and N\phi. More of a sketch than a proof, let me know if you would like me to try and present it clearer!