Consider this model that could generate correlated Poisson variables. Let $Y$, $Y_1$ and $Y_2$ be three independent Poisson variable with parameters $r$, $\lambda_1$ and $\lambda_2$. Let $X_i=Y_i+Y$ for $i=1,2$. Then $X_1$ and $X_2$ are both Poisson with parameters $\lambda_1$ and $\lambda_2$. They have the correlation $\rho=\frac{r}{\sqrt{(\lambda_1+r)(\lambda_2+r)}}$ Now the joint distribution can be derived as $P[X_1=i,X_2=j]=e^{-(r+\lambda_1+\lambda_2)}\sum_{k=0}^{i\wedge j}\frac{r^k}{k!}\frac{\lambda_1^{(i-k)}}{(i-k)!}\frac{\lambda_2^{(j-k)}}{(j-k)!}$ The case for a bivariate Poisson process is immediate from here.
You could look at the Johnson and Kotz book on multivariate discrete distributions for more information (this construction of a bivariate Poisson distribution is not unique). Also, it has the drawback that $\rho \in [0, \min(\lambda_1, \lambda_2)/\sqrt{\lambda_1\lambda_2} ]$ when $\lambda_1 \neq \lambda_2$ as discussed by Genest et al. 2018.