How is the Mann Whitney U test statistic derived? The t-test statistic is derived from the CLT and allows for some probabilistic statements to be made.
$$ U_1 = n_1(n_1+1)/2 + n_1 n_2 - R_1 .$$
I know the test does not assume a probability distribution; though, I was wondering where this formula came from, why it is defined in this way and how we are able to make probability statements about the similarity of two groups from it.

A general outline is worth to delve into:
The Mann-Whitney U statistic can be defined (cf. $\rm [I]$) as the number of times $Y$ exceeds $X$ in the combined ordered arrangement of independent samples $\langle X_i\rangle_{i=1}^m, ~\langle Y_j\rangle_{j=1}^n.$ Succinctly
$$U:= \sum_{i=1}^m\sum_{j=1}^n\boldsymbol 1_{X_i< Y_j}\tag 1\label 1.$$
Consider now the Wilcoxon sum statistic
$$W:= \sum_{j=1}^n R(Y_j),\tag 2\label 2$$
where $R(Y_j)$ is the rank of $Y_j$ in the combined sample. How are $\eqref 1, ~\eqref 2$ related (cf. $\rm [II]$)?
For that, notice
$$R(Y_j) = \sum_{i=1}^m \boldsymbol 1_{X_i< Y_j} + \sum_{k=1}^n\boldsymbol 1_{Y_k< Y_j}+1;\tag 3\label 3$$ substituting $\eqref 3$ in $\eqref 2$ yields
\begin{align}W &= \left(\sum_{j=1}^n\sum_{i=1}^m\boldsymbol 1_{X_i< Y_j}\right) + \left(\sum_{j=1}^n\sum_{k=1}^n\boldsymbol 1_{Y_k< Y_j}+n\right)\\&= U +\frac{n(n+1)}{2}\tag 4\label 4.
\end{align}
So, from $\eqref 4,$ the Mann-Whitney U statistic can be equivalently expressed as $$U = W-\frac{n(n+1)}{2}.\tag 5\label 5$$

In adherence to the notations used by OP, consider the two samples of sizes $n_1, ~n_2$ and the sum of the corresponding ranks of the observations in the respective samples be denoted by $R_1, ~R_2.$ Also, let $N:= n_1+n_2.$
Naturally,
$$R_1 + R_2 = \frac{N(N+1)}2. $$ That is,
\begin{align}R_1+R_2&= \frac{(n_1+n_2)(n_1+n_2+1)}2\\ &= \frac{n_1^2+ n_1n_2+n_1 + n_2n_1+ n_2^2 +n_2}{2}\\&= \frac{n_1(n_1+1)}{2} + n_1n_2 + \frac{n_2(n_2+1)}{2}\tag{6.a}\label{a}; \end{align}
from $\eqref a,$
$$R_1 - \frac{n_1(n_1+1)}{2} =n_1n_2 + \frac{n_2(n_2+1)}{2} - R_2;\tag{6.b}\label b  $$
in $\eqref b,$ LHS is of the form $\eqref 5;$ hence
$$U_2 = n_1n_2 + \frac{n_2(n_2+1)}{2} - R_2.$$
Similarly,
$$U_1 = n_1n_2 + \frac{n_1(n_1+1)}{2} - R_1.$$ Interestingly, it follows $U_1 + U_2 = n_1n_2.$

References:
$\rm[I]$ Nonparametric Statistical Inference, Jean Dean Gibbons, Subhabrata Chakraborti, Marcel Dekker, $2003,$ sec. $6.6,$ p. $268.$
$\rm [II]$ Nonparametric Statistical Methods, Myles Hollander, Douglas A. Wolfe, Eric Chicken,  John Wiley & Sons, $2014,$ sec. $4.1,$ pp. $126-127.$
$\rm [III]$ Biostatistical Analysis, Jerrold H. Zar, Pearson Education, $2014,$ p. $171.$
