So after further work on the problem (and crucially asking a PhD for orientation), I believe I have solved the problem.
Part 1: General method using the example of $1^1\Sigma_u$ -> $1^1\Sigma_g$
$<1^1\Sigma_g|q_1r_1+q_2r_2|1^1\Sigma_u>$ = $\int{[A(1)B(2)+A(2)B(1)](q_1r_1+q_2r_2)[A(1)A(2)-B(1)B(2)]}dr_1dr_2$
I will leave out the q's from now on, except where the final result is written. The normalization coefficients are not necessary to include as will be shown later.
For this $1^1\Sigma_u\>$ -> $1^1\Sigma_g$ transition, I had the following 8 terms after multiplying out the integrand:
+$\int{[r_1A(1)^2A(2)B(2)]dr_1dr_2}$
+$\int{[r_1A(2)^2A(1)B(1)]dr_1dr_2}$
+$\int{[r_2A(1)^2A(2)B(2)]dr_1dr_2}$
+$\int{[r_2A(2)^2A(1)B(1)]dr_1dr_2}$
-$\int{[r_1A(1)B(1)B(2)^2]dr_1dr_2}$
-$\int{[r_1A(2)B(2)B(1)^2]dr_1dr_2}$
-$\int{[r_2A(1)B(1)B(2)^2]dr_1dr_2}$
-$\int{[r_2A(2)B(2)B(1)^2]dr_1dr_2}$
These can be simplified by recognizing that a non-diverging, normalized function integrates to 1, eg. $\frac{1}{\sqrt{\pi}}\int{e^{-x^2}dx}=1$. Of course I left out the normalization coefficients of the eigenstates above (as well as below here) but adding them in should by definition yield that result.
Thus 4 terms can be immediately reduced:
+$\int{[r_1A(1)^2A(2)B(2)]dr_1dr_2}$
+$\int{[r_1A(2)^2A(1)B(1)]dr_1dr_2}$ = $\int{A(2)^2dr_2}\int{[r_1A(1)B(1)]dr_1}$ = $1\int{[r_1A(1)B(1)]dr_1}$
+$\int{[r_2A(1)^2A(2)B(2)]dr_1dr_2}$ = $\int{A(1)^2dr_1}\int{[r_2A(2)B(2)]dr_2}$ = $1\int{[r_2A(2)B(2)]dr_2}$
+$\int{[r_2A(2)^2A(1)B(1)]dr_1dr_2}$
-$\int{[r_1A(1)B(1)B(2)^2]dr_1dr_2}$ = -$\int{B(2)^2dr_2}\int{[r_1A(1)B(1)]dr_1}$ = -$1\int{[r_1A(1)B(1)]dr_1}$
-$\int{[r_1A(2)B(2)B(1)^2]dr_1dr_2}$
-$\int{[r_2A(1)B(1)B(2)^2]dr_1dr_2}$
-$\int{[r_2A(2)B(2)B(1)^2]dr_1dr_2}$ = -$\int{B(1)^2dr_1}\int{[r_2A(2)B(2)]dr_2}$ = -$1\int{[r_2A(2)B(2)]dr_2}$
Clearly these four simplified terms nicely cancel out, leaving us with the four remaining terms. At this point, we define the overlap integral $S_{AB}=\int{A(x)B(x)dx}$ to reduce the remaining terms.
+$\int{[r_1A(1)^2A(2)B(2)]dr_1dr_2}$ = $S(2)_{AB}\int{r_1A(1)^2dr_1}$
+$\int{[r_2A(2)^2A(1)B(1)]dr_1dr_2}$ = $S(1)_{AB}\int{r_2A(2)^2dr_2}$
-$\int{[r_1A(2)B(2)B(1)^2]dr_1dr_2}$ = -$S(2)_{AB}\int{r_1B(1)^2dr_1}$
-$\int{[r_2A(1)B(1)B(2)^2]dr_1dr_2}$ = -$S(1)_{AB}\int{r_2B(2)^2dr_2}$
To proceed further, we need to establish a relation for $\int{r_1A(1)^2dr_1}$ and $\int{r_1B(1)^2dr_1}$. For this I considered the symmetry of the system:
You can prove such an integral is zero by setting up the equation as:
$\int{xe^{-(x+c)^2}e^{-(x-c)^2}dx}$
or
$\int{(x+\frac{c_1+c_2}{2})e^{-(x+c_1)^2}e^{-(x+c_2)^2}dx}$
for a system not at the origin.
Steps for solving the latter integral
You multiply out the exponents, pull out the constants, then complete the square in the exponent of the exponential function by implicitly adding and subtracting $\frac{c_1+c_2}{4}$. Pulling out the leftover constant from the implicit subtraction, you can substitute in $y=\sqrt{b}(x+\frac{c_1+c_2}{2})$ for $e^{-b(x+\frac{c_1+c_2}{2})^2}$. From here, integration by parts with $v$ = $\int{e^{-y^2}dy}=\sqrt{\pi}$ will lead to 0.
Because we have now shown this is an odd function, we know that the integral to the left of the origin in the plot above must be equal and negative to the integral on the right. This means that
$\int{r_1A(1)^2dr_1}=-\int{r_1B(1)^2dr_1}$
and substituting this in for our 4 terms above will allow us to add together $r_1$ and $r_2$ terms, respectively. We are now left with just two terms from our original problem:
$2S(2)_{AB}\int{r_1A(1)^2dr_1}$
$2S(1)_{AB}\int{r_2A(2)^2dr_2}$
and this is in fact our solution. In summary:
$<1^1\Sigma_g|q_1r_1+q_2r_2|1^1\Sigma_u>$ = $q_1S(2)_{AB}\int{r_1A(1)^2dr_1}+q_2S(1)_{AB}\int{r_2A(2)^2dr_2}$
where the normalization constants of $\frac{1}{\sqrt{2}}$ even do us the favor of canceling out the 2 coefficient. Our transition dipole moment is non zero, as expected. At large distances however, it will disappear as the overlap integral goes to zero.
Part 2: $2^1\Sigma_g$ -> $1^1\Sigma_g$
$<1^1\Sigma_g|q_1r_1+q_2r_2|2^1\Sigma_g>$ = $\int{[A(1)B(2)+A(2)B(1)](q_1r_1+q_2r_2)[A(1)A(2)+B(1)B(2)]}dr_1dr_2$
For the ungerade -> gerade transition above, we started with 8 terms after multiplying out the integrand, 4 of which were negative. In this gerade -> gerade transition, all of those terms will now be positive.
We began our solution above by removing normalized integrals, thereby canceling out 4 terms. This will no longer be possible for the gerade-gerade transition, since there are no negative terms to subract. However, we later showed that we are dealing with an odd function. By this symmetry, these same 4 integrals will still go to zero. Subtraction is not necessary to cancel them out.
The remaining 4 terms from our ungerade -> gerade transitoin were added together by the relationship
$\int{r_1A(1)^2dr_1}=-\int{r_1B(1)^2dr_1}$
In the case of a gerade-gerade transition, there is no negative sign to cancel out this antisymmetry. This means the last 4 terms will not add together but subtract, leading to zero.
Thus all 8 terms result in zero transition dipole moment, as expected of a gerade-gerade transition.
Part 3: $2^1\Sigma_g$ -> $1^1\Sigma_u$
Finally, we have
$<1^1\Sigma_u|q_1r_1+q_2r_2|2^1\Sigma_g>$ = $\int{[A(1)A(2)-B(1)B(2)](q_1r_1+q_2r_2)[A(1)A(2)+B(1)B(2)]}dr_1dr_2$
The result after multiplying out was
$\int(r_1+r_2)[A(1)^2A(2)^2-B(1)^2B(2)^2]dr_1dr_2$
This can be broken down into 2 terms in the same way we did above. First, we utilize their normalization to remove any integrals of the form $\int{C(x)^2dx}$ as 1. Then we combine A and B terms by their symmetry $\int{A(x)dx}$ = -$\int{B(x)dx}$. These steps leave us with our two terms:
$\int{r_1A(1)^2dr_1}$
$\int{r_2A(2)^2dr_2}$
where the 2 coefficient has also again been canceled out by the normalization constant on the eigenstates. This transition is clearly non-zero as expected of a gerade to ungerade transition. The solution,
$<1^1\Sigma_u|q_1r_1+q_2r_2|2^1\Sigma_g>$ = $q_1\int{r_1A(1)^2dr_1}+q_2\int{r_2A(2)^2dr_2}$
has no overlap integral. Therefore, the transition dipole moment increases linearly with respect to separation, at a slope of 2 (for $q_1=q_2=2e$ in this doubly ionized species). This confirms a result I was getting with CI calculations of the transition dipole moment via GameSS.
Conclusions
- Solved the transitions and reached the expected results. The transition moment is non-zero for gerade -> ungerade and zero for gerade -> gerade. An overlap integral as well as symmetry between the atoms were key points in the solution.
- The overlap integral for the $1^1\Sigma_u$ -> $1^1\Sigma_g$ transition means that at long distances, $S_{AB}$ goes to zero and the dipole moment for this transition disappears; transitions at long distances are not possible between these eigenstates.
- The lack of an overlap integral for the $2^1\Sigma_g$ -> $1^1\Sigma_u$ transition means that the dipole moment is proportional to the distance, and even at large distances it does not disappear but increases. This may have to do with their being energetically asymptotic with respect to separation, as they both represent the $He^{2+}+He(1s^2)$ configuration.
I hope that helped someone!