Suppose that for two discrete random variables $X_1$ and $X_2$, we know their conditional distributions. Namely $$X_1~|~X_2 = x_2 \sim \mathrm{Poisson}(\lambda_1 + ax_2),$$ $$X_2~|~X_1 = x_1 \sim \mathrm{Poisson}(\lambda_2 + bx_1).$$

We want to calculate their joint distribution, $p(X_1, X_2)$.

My own idea is to divide the equations above and find $\frac{p_{X_1}(x_1)}{p_{X_2}(x_2)}$ and then sum over all values of $X_1$ to find $1/p_{X_2}(x_2)$. But for this, I have to compute a very bad series.

Do you have any idea for this problem?

P.S. The series I have to compute is of form $$\sum_{n=0}^\infty \frac{c^n}{n!(a+nb)^k},$$ which I hardly believe to have a closed form.