In a sequence of Bernoulli trials with success chance $p$, the probability that we will see $k$ failures before the $r$-th success is given by the negative binomial distribution (NBD) with pmf: $$ Pr(k) = \binom{k+r-1}{r-1} p^r (1-p)^k $$ But what if we had two sequences of simultaneous Bernoulli trials and we wanted to know the probability of how many trials we need, before we see at least $r_1$ successes in the first sequence and $r_2$ successes in the second, with success chances $p_1$ and $p_2$, respectively?
For simplicity, consider the case where $r_1 = r_2 = 1$. The NBD pmf simplifies to: $$ Pr(k) = p (1-p)^k $$ Let's distinguish two independent random variables $X$ and $Y$ with their respective pmf's: $$ Pr_x(k_x) = p_x (1-p_x)^{k_x} \quad \text{and} \quad Pr_y(k_y) = p_y(1-p_y)^{k_y}$$
If we perform a sequence of trials for both $X$ and $Y$ simultaneously, and it takes $n$ trials before we see at least one success in each sequence, then at least one sequence must have had its first success exactly after $n$ trials (or $n-1$ failures). Suppose that this was the $X$-sequence. Then the $Y$-sequence could also have had its first success after $n-1$ failures. The probability of this happening is: $$ Pr_x(n-1) \cdot Pr_y(n-1) $$ But the $Y$-sequence could also have seen its first success before the $n$-th trial. The probability of that happening is: $$ Pr_x(n-1) \cdot Pr_y(< n-1) = Pr_x(n-1) \cdot \Bigg(\sum_{j = 0}^{n-2} Pr_y(j)\Bigg) $$ Lastly, there is the case where the $Y$-sequence has its first success after exactly $n-1$ failures, while the $X$-sequence had its first success before the $n$-th trial. The probability of this occurring is: $$ Pr_x(< n-1) \cdot Pr_y(n-1) = \Bigg(\sum_{j = 0}^{n-2} Pr_x(j)\Bigg) \cdot Pr_y(n-1) $$ The summations in the two previous expressions are geometric. As long as we assume that $p_x, p_y \neq 0$, we can simplify them. Let $(q_x, q_y) := (1-p_x, 1-p_y)$ and by adding all three probabilities together we find a nice, symmetric expression: $$ Pr(n) = p_x q_x^{n-1} p_y q_y^{n-1} + p_x q_x^{n-1} (1 - q_y^{n-1}) + p_y q_y^{n-1} (1 - q_x^{n-1}) $$ I've verified this expression with a few experiments in Python. Is the constructive argument I've given above enough to prove that the pmf is correct?
So, what do we do in the more general case where we're looking for more than one success? I've already struggled quite a bit with this simple case. And more importantly, is this sort of distribution well-known and if so, what is it called?