Why is $ab$ likely to have more divisors than $(a-b)(a+b)$? Consider the two numbers $ab$ and $(a-b)(a+b), \gcd(a,b) = 1, 1 \le b < a$. On an average, which of these two numbers has more distinct prime factors? All the prime factors of $a$ and $b$ divide $ab$ and similarly all the prime factors of $a-b$ and $a+b$ divide $(a-b)(a+b)$. So one number does not seem to have an obvious advantage over the other. However if we look at the data than we see that $ab$ dominates.
Let $f(x)$ be the average number of distinct prime factors in all such $ab, a \le x$ and
$g(x)$ be the average number of distinct prime factors in all such $(a-b)(a+b), a \le x$.
Update: Experimental data for the first $6.1 \times 10^{9}$ pairs of $(a,b)$ shows that $f(x) - g(x) \sim 0.30318$. Instead of distinct prime factors, if we count the number of divisors than $f(x) - g(x) \sim 0.848$.
Question: Why is $ab$ likely to have more distinct prime factors or divisors than $(a-b)(a+b)$ and what is the limiting value of $f(x) - g(x)$?
 A: Note if $\gcd(a, b) = 1$, then $\gcd(a, a - b) = 1$. Also, $\gcd(a - b, a + b) = d$ where $d = 1$ or $d = 2$.
Fix any $a \gt 2$, and let $m = \left\lceil \frac{a}{2} \right\rceil - 1$. Consider each $1 \le b_1 \le m$ where $\gcd(a, b_1) = 1$. We have $b = b_1$ giving $a(b_1)$ and $(a - b_1)(a + b_1)$. Correspondingly, $b = a - b_1$ gives $a(a - b_1)$ and $b_1(2a - b_1)$. Note over the specified range of $b_1$, these $2$ groups represents the values being checked for every coprime $b$ exactly once.
For simpler algebra, let $h(x)$ be the # of distinct odd prime factors of $x$. Thus, the total # of odd prime factors of $ab$ over the coprime $b$ is
$$\begin{equation}\begin{aligned}
q(a) & = \sum_{b_1=1,\,\gcd(a,b_1)=1}^{m}(h(a(b_1)) + h(a(a - b_1))) \\
& = \sum_{b_1=1,\,\gcd(a,b_1)=1}^{m}(2h(a) + h(b_1) + h(a - b_1))
\end{aligned}\end{equation}\tag{1}\label{eq1A}$$
Similarly, the total # of odd prime factors of $(a - b)(a + b)$ over the coprime $b$ is
$$\begin{equation}\begin{aligned}
r(a) & = \sum_{b_1=1,\,\gcd(a,b_1)=1}^{m}(h((a - b_1)(a + b_1)) + h(b_1(2a - b_1)) \\
& = \sum_{b_1=1,\,\gcd(a,b_1)=1}^{m}(h(a - b_1) + h(a + b_1)) + h(b_1) + h(2a - b_1))
\end{aligned}\end{equation}\tag{2}\label{eq2A}$$
Next, we have
$$s(a) = q(a) - r(a) = \sum_{b_1=1,\,\gcd(a,b_1)=1}^{m}(2h(a) - h(a + b_1) - h(2a - b_1)) \tag{3}\label{eq3A}$$
Note $a + b_1 \gt a$ and $2a - b_1 \gt a$. If $a$ is odd and $b_1$ is even then $a + b_1$ and $2a - b_1$ are odd and, on average, would have somewhat more distinct odd prime factors than $a$. However, with $a$ odd and $b_1$ odd, then $a + b_1$ and $2a - b_1$ would both be even, so their odd parts (i.e., largest odd factors) would be smaller than $a$, so they would likely have fewer distinct odd prime factors than $a$, with this at least roughly canceling the effect of when $b_1$ is even to make $s(a)$ overall close to $0$. Finally, if $a$ is even, then $b_1$, $a + b_1$ and $2a - b_1$ are all odd, so since the odd part of $a$ is considerably less than $a$, the $s(a)$ would generally be negative on average. Although the number of cases of $b_1$ here would only be somewhat more than half of the ones for $a$ being odd (see the next paragraph for details), the total effect should be for $s(a)$ to be negative when summed over many $a$.
Regarding the occurrences of the prime factor of $2$, for even $a$, let $a = 2^{\nu_2(a)}c$ where $c$ is odd. Then by the multiplicative property of Euler's totient function, $e = \frac{\varphi(a)}{a} = \frac{2^{\nu_2(a) - 1}\varphi(c)}{2^{\nu_2(a)}c} = \left(\frac{1}{2}\right)\left(\frac{\varphi(c)}{c}\right)$. In comparison, with $f = \frac{\varphi(a+1)}{a + 1}$, since $a + 1 \gt c$ then, on average, $a + 1$ would have more distinct prime factors and, thus, $f \lt \frac{\varphi(c)}{c} \; \to \; e \gt \frac{f}{2}$. Therefore, the number of $b$ when $a$ is even on average would be somewhat more than half of that of $a + 1$ (unfortunately, I'm not sure how to calculate how much more). Thus, the case of $a$ being odd with $b$ being odd and $b$ being even are equal (since $b_1$ and $a - b_1$ have opposite parities), but the case of $a$ being even and $b$ being odd is somewhat more common.
For the second & third cases (i.e., $a$ & $b$ having different parities), $ab$ is even (so there's somewhat more than a $\frac{2}{3}$ chance), while $(a - b)(a + b)$ is even only for the first case of $a$ & $b$ both being odd (i.e., a somewhat less than a $\frac{1}{3}$ chance), giving a difference of somewhat more than $\frac{1}{3}$. Accounting for the result in \eqref{eq3A} being generally a small negative, this means that $f(x) - g(x)$ would be somewhat smaller (although I'm not sure how to accurately calculate the amount less), and I believe it'll make the overall result somewhat less than $\frac{1}{3}$. Thus, I suspect your latest result (from your comment) of $f(x) - g(x) \sim 0.3016$ for the first $1.4 \times 10^{9}$ pairs of $(a,b)$ is likely relatively close to the asymptotic value.
A: On one hand,
\begin{align*}
\sum_{a,b\le x} \omega(ab) &= \sum_{a,b\le x} \sum_{\substack{p\le x \\ p\mid ab}} 1 = \sum_{p\le x} \sum_{\substack{a,b\le x \\ p\mid ab}} 1 \\
&= \sum_{p\le x} \biggl( \sum_{\substack{a,b\le x \\ p\mid a}} 1 + \sum_{\substack{a,b\le x \\ p\mid b}} 1 - \sum_{\substack{a,b\le x \\ p\mid a,\, p\mid b}} 1 \biggr) \\
&= \sum_{p\le x} \biggl( \bigl( \tfrac xp+O(1) \bigr)(x+O(1)) + (x+O(1))\bigl( \tfrac xp+O(1) \bigr) - \bigl( \tfrac xp+O(1) \bigr)^2 \biggr) \\
&= 2x^2 \sum_{p\le x} \tfrac1p - x^2 \sum_{p\le x} \tfrac1{p^2} + O\biggl( x \sum_{p\le x} \bigl( 1+\tfrac1p \bigr) \bigg) = 2x^2 \sum_{p\le x} \tfrac1p - x^2 \sum_p \tfrac1{p^2} + o(x^2).
\end{align*}
On the other hand, $p\mid(a+b)$ and $p\mid(a-b)$ simultaneously if and only if either $p\mid a$ and $p\mid b$, or $p=2$ and $a$ and $b$ are both odd. Therefore
\begin{align*}
\sum_{a,b\le x} \omega\bigl( (a+b)(a-b) \bigr) &= \sum_{a,b\le x} \sum_{\substack{p\le x \\ p\mid (a+b)(a-b)}} 1 = \sum_{p\le x} \sum_{\substack{a,b\le x \\ p\mid (a+b)(a-b)}} 1 \\
&= \sum_{p\le x} \biggl( \sum_{\substack{a,b\le x \\ p\mid (a+b)}} 1 + \sum_{\substack{a,b\le x \\ p\mid (a-b)}} 1 - \sum_{\substack{a,b\le x \\ p\mid (a+b),\, p\mid (a-b)}} 1 \biggr) \\
&= \sum_{p\le x} \biggl( \sum_{a\le x} \sum_{\substack{b\le x \\ b\equiv-a\,(\mathop{\rm mod}\,p)}} 1 + \sum_{a\le x} \sum_{\substack{b\le x \\ b\equiv a\,(\mathop{\rm mod}\,p)}} 1 - \sum_{\substack{a,b\le x \\ p\mid a,\, p\mid b}} 1 \biggr) \\
&\qquad{}- \sum_{\substack{a,b\le x \\ a,b \text{ both odd}}} 1 \\
&= \sum_{p\le x} \biggl( (x+O(1))\bigl( \tfrac xp+O(1) \bigr) + (x+O(1))\bigl( \tfrac xp+O(1) \bigr) - \bigl( \tfrac xp+O(1) \bigr)^2 \biggr) \\
&\qquad{}- \bigl( \tfrac x2+O(1) \bigr)^2 \\
&= 2x^2 \sum_{p\le x} \tfrac1p - x^2 \sum_p \tfrac1{p^2} - \tfrac{x^2}4 + o(x^2).
\end{align*}
From these two asymptotic formulas it follows that the difference of the two sums is asymptotic to $\frac{x^2}4$, so that the average difference in the number of distinct prime factors is asymptotically $\frac14$.
Heuristically (in hindsight), the difference is entirely caused by the prime $2$: since $a$ and $b$ are even or odd independently, there is a $\frac34$ chance that $p=2$ will contribute to $\omega(ab)$; but since $a+b$ and $a-b$ are both even or both odd, there is only a $\frac12$ chance that $p=2$ will contribute to $\omega\bigl( (a+b)(a-b) \bigr)$.
A: In case A, looking at the set $ \{(1,0),(0,1),(1,1)\} $, the probability that 2|ab is 2/3
In case B, looking at the same set, the probability that $ 2|(a-b)(a+b) $ is $ 1/3
 $
For all other odd primes, the probability that $p|ab$ is the same as $p|(a-b)(a+b)$, and is: $2(p-1)/(p^2-1) = 2/(p+1)$
However, for large primes, $a<p<2a$, it is possible that $p|(a+b)$ but never $p|ab$.
Hence, the result that comes out empirically is somewhat less than 1/3.
At 1/3 of the average, some of the a+b=p mentioned must be subtracted.
