Consider the function $F: (\lambda_1,\lambda_2) \rightarrow (a,c)$, where $a = P(X_1 > X_2)$ and $c = 1 - a - b = P(X_1 < X_2)$. We claim that $F$ is invertible, considered as a map from $(0,\infty)^2$ to the open triangular set $T = \{(x,y) \in \mathbb{R}^2; x>0, y>0, x + y < 1\}$.

**$F$ is injective**

First, we compute the Jacobian.
\begin{eqnarray}
\frac{da}{d\lambda_1} &=&-a + e^{-\lambda_1-\lambda_2}\sum_{k = 0}^\infty\sum_{j > k}\frac{\lambda_1^{j-1}\lambda_2^k}{(j-1)!k!}\\&=&-a + e^{-\lambda_1-\lambda_2}\sum_{k = 0}^\infty\sum_{j \geq k}\frac{\lambda_1^j\lambda_2^k}{j!k!}\\
&=&-P(X_1>X_2) + P(X_1\geq X_2)\\
&=&P(X_1 = X_2).
\end{eqnarray}
Similar calculations for the other partial derivatives gives:
$$
J = \left(\begin{array}{cc}P(Z = 0) &-P(Z = 1)\\-P(Z = -1) &P(Z=0)\end{array}\right),
$$
where $Z = X_1 - X_2$ is the Skellam distributed random variable. The probability mass function for $Z$ is given by $P(Z=k) = e^{-\lambda_1 -\lambda_2}(\lambda_1/\lambda_2)^{k/2}I_{|k|}(2\sqrt{\lambda_1\lambda_2})$, where $I_{k}$ is the modified Bessel function of the first kind.

Now, $J$ is a P-matrix, i.e. all its principal minors are positive. Indeed, $P(Z=0) > 0$, and the determinant $P(Z=0)^2 - P(Z=1)P(Z=-1)$ equals $e^{-2\lambda_1 -2\lambda_2}(I_0(2\sqrt{\lambda_1\lambda_2})^2-I_1(2\sqrt{\lambda_1\lambda_2})^2)$, which is positive since $I_0 > I_1$.

But now we have everything we need to apply the "Fundamental Global Univalence Theorem (Gale-Nikaido-Inada)" (see On Global Univalence Theorems, Lecture Notes in Mathematics Volume 977, 1983, pp 17-27), which gives us injectivity of $F$.

**$F$ is surjective**

It is clear by construction that $F = (a,c)$ does not take values outside the closure of $T$: $a$ and $c$ are probabilities of disjoint events, so they are non-negative and sum to at most $1$. Moreover, the range of $F$ is simply connected, since $F$ is continuous and its domain is simply connected. Consider the three line segments that constitute the boundary of $T$: $R_1 = \{(x,0); 0 < x < 1\}$, $R_2 = \{(0,y); 0 < y < 1\}$, and $R_3 = \{(x,y); x>0,y>0,x + y = 1\}$. We are done if we can show that $F$ attains values arbitrarily close to each point in $R_1\cup R_2\cup R_3$, but never takes values in that set. (One should also consider the three corners, $(0,0)$, $(1,0)$ and $(1,0)$, but the arguments are similar for these cases.)

Consider a point $(x,0) \in R_1$. It is clear that this value is never attained by $F$, since $P(X_2 > X_1) = 0$ would imply either vanishing $\lambda_2$, or infinite $\lambda_1$. However, values arbitrarily close to $(x,0)$ are attained: choose $\lambda_1$ such that $P(X_1 > 0) = x$. Then $F(\lambda_1,\lambda_2) \rightarrow (x,0)$ as $\lambda_2 \rightarrow 0$. The set $R_2$ is treated similarly.

Finally, consider a point $(x,y) \in R_3$. By symmetry, we can assume that $x > y$. This point can not be in the range of $F$, since that would imply $P(X_1 = X_2) = 0$, which is impossible for finite $\lambda_1$ and $\lambda_2$. Now, $1 - a - c \rightarrow 0$ when $\lambda_2 \leq \lambda_1$ and $\lambda_2\rightarrow \infty$, so we can choose $\lambda_2$ sufficiently large and send $\lambda_1$ from $\lambda_1 = \lambda_2$ towards infinity: $F(\lambda_1,\lambda_2)$ will describe a curve from a point in the epsilon ball $B_\varepsilon(0.5,0.5)$ to a point in $B_\varepsilon(1,0)$, always staying at most $\varepsilon$ away from $R_3$, hence passing within $\varepsilon$ of $(x,y)$.