Expected size of least prime factor of $(16^p-1)/15$ Let's define a function $g(x)$ as follows:
$g(x)=-x+\displaystyle\sum_{\substack{3\le p\le x/2}\\{p \text{ prime}}}{\log_2 \left( \text{LeastPrimeFactor} \left( \frac{16^p-1}{15} \right) \right) }$
This function comes from a specific rabbit-hole of investigation around the Freudenthal Problem.  I want to know if $\liminf g(x) > 0$.  This leads directly to the question in the title.
The sum ranges over odd primes less than $x/2$, so the crux of the problem is finding the factors of Mersenne numbers, which makes me immediately think this is out of reach.  However, I've seen some heuristics discussed that make it seem like there might be a way forward.
Naively, we know 4 factors of $16^p-1$:

*

*$2^p-1$

*$2^p+1$

*$2^p-2^\frac{1+p}{2}+1$

*$2^p+2^\frac{1+p}{2}+1$
These are all of order $2^p$, one is divisible by 5 and another by 3, but there are not guaranteed factors beyond this.  Plugging them in gives us the sum of all odd primes up to $x/2$, dominating the $-x$ and leading to $+\infty$.
On the other hand, for $p\ge 7$ we know the lower bound for the least prime factor of this expression is $2p+1$, a bound which is attained quite often.  Summing $\log_2 \left( 2p+1 \right)$ gives something $O(x)$, and if we replace it throughout, $g(x)$ tends to $-\infty$.
Using actual data, mostly scraped from Will Edgington's database and assembled in Mathematica, it looks like the big factors eventually start winning.  Of course, I'm hesitant to make claims using only data, particularly when discussing not-yet-found factors.  For example, when $p=2267$, it seems that we know one prime factor with 40 digits and haven't proven that it is the least prime factor.  But when $p=11423$, we don't know any prime factors!

Is there some reasonable way to predict in-bulk, the expected behavior of the least prime factors of these Mersenne-adjacent numbers?  If so, what does it predict about $g(x)$?
Edit:
I tried my hand at a probabilistic number theory approach.
The progression $2np+1$ includes all factors of $(16^p-1)/15$.
The size of the $k$th prime in that progression goes as $2pk \log(2pk) / \log(2p)$, but the log of that value can be approximated reasonably as $\log(2pk)$.
The upper limit of $k$ is $2^p$ due to the known factors.  You could also put $2^{p/2}$ here and the expressions below would not change at leading order.
Some integrals show that

*

*$\int^{2^p} \frac{1}{k} \,dk \rightarrow p \log 2$

*$\int^{2^p} \frac{1}{k \log k} \,dk \rightarrow \log p$

*$\int^{2^p} \frac{1}{k \log^2 k} \,dk \rightarrow 1/\log 2$.

If I follow the argument in Expected smallest prime factor with appropriate generality, I think the probability the $k$th prime in the progression is the smallest prime factor, goes as $\frac{1}{k \log^2 k}$.
This gives rise to an expectation value for $\log_2 \left( \text{LeastPrimeFactor} \left( 2^p-1 \right) \right)$ that is $(1 + 1/\log2)\log p + O(\log \log p)$, and a variance that is $O(p)$.
This might be a classic result?  Or it might be an erroneous result.  But converting the base it sounds like the expected number-of-bits in the smallest factor of the $p$th Mersenne number is $1 + \log 2 = 1.69315..$ times the number of bits in $p$ itelf.
The variance is large, and the distribution has a hard minimum of $\log(2p+1)$, so among many samples of this distribution there will be cases where the log of the factor is itself $O(\sqrt{p})$.
Now, the log-least-prime-factor of $(16^p - 1)/15$ would be the minimum of 4 reasonably-independent draws from this distribution, so it avoids the $O(\sqrt{p})$ draws with higher likelihood.  But by this point I feel like the approximations and assumptions are going to catch up with me before I get anywhere.
So, I'd like to have help with continuing/refining this argument until it can predict something reasonably useful about $g(x)$.
 A: Ok, I have at least some progress worth sharing.
Put $f(p) = \log_2\left(\text{LeastPrimeFactor}\left(\left(16^p-1\right)/15\right)\right)$.  Let's treat $f(p)$ as a random variable, and compute its expectation value as a function of $p$:
To do so, consider "some value of order $2^p$ which is divisible only by numbers of the form $2pn+1$".  Namely, this definition encompasses Mersenne numbers of prime exponent $2^p-1$, as well as $(2^p+1)/3$, $2^p\pm 2^{(1+p)/2}+1$, and $(2^p\mp 2^{(1+p)/2}+1)/5$.
The product of these 4 expressions is our current target number, and because $2^p+1-(2^p-1)=2$, $4^p+1-(4^p-1)=2$, and $2^p + 2^{(1+p)/2}+1 - (2^p - 2^{(1+p)/2}+1) = 2^{(3+p)/2}$, there is no way for any pair of them to have factors in common.
Among the progression $2pn+1$, which numbers can be the smallest prime factor?  Only the primes, for one.  Additionally, they must be either less than $2^{p/2}$ (otherwise the cofactor would be smaller and thus have smaller prime factors), or be the expression itself.  With our goal here, we will be able to neglect the case where all 4 expressions are prime.  It will make $f(p)$ larger than we predict, if it does happen for some large $p$, but we are able to make progress whether or not it ever happens for $p>7$.
By the behavior of PNT in arithmetic progressions, we expect the $k$th prime $p_k$ in the progression $2pn+1$ to be approximately $p_k \approx pk\log(2pk)$.  This holds up well vs data.
The probability $\pi_k$ that $p_k$ is the smallest prime factor of $(16^p-1)/15$ is the product of:

*

*The probability that it divides one of the 4 known factors of order $2^p$. If we just used sieve theory, this would be $4/p_k$ for sufficiently large $p$.  However, we will add a factor $A$ in as well, for reasons understood later.

*The probability that for all $k' < k$, The $k'$th prime is not a divisor of any of the 4 known factors. From the above, this is a product $\prod_{k'<k}{\left(1-4A/p_{k'}\right)}$.

When $p$ and $k$ are sufficiently large,
$\frac{\pi_{k+1}}{\pi_{k}} = \frac{k\log(2pk)}{(k+1)\log(2pk+2p)}\left(1-\frac
{4A}{pk\log(2pk)}\right) \approx 1-\frac{4A+p+p\log(2pk))}{pk\log(2pk)} \rightarrow \partial_k(\log(\pi_k))=-\frac{4A+p+p\log(2pk))}{pk\log(2pk)}$
Which suggsts that
$\pi_k \propto \exp\left(-\int\frac{4A+p+p\log(2pk))}{pk\log(2pk)}\right) \approx k^{-1}\log(2pk)^{-(1+4A/p)}$
For any $A>0$ this is a convergent sum, but most of its support at large $p$ is over $k>2^p$, which is too large to be a factor.  We need an upper bound on $k$ at $k_{\text{max}}=p^{-2}\times2^{p/2}$, so that $p_{k_\text{max}} \approx 2^{p/2}$.  We cannot simply truncate the sum and then divide throughout by the partial sum, as this changes the "probability that $p_k$ divides the expression" without corresponding changes to "the probability that for all $k'<k$, $p_k'$ does not divide the expression".  This is why we have the factor $A$.  We pick $A$ such that the partial sum is self-consistent, as follows:
$\pi_k \approx \frac{4A}{pk \log(2pk)}\left(\frac{\log(2p)}{\log(2pk)}\right)^{4A/p}\hspace{15pt}\sum_{k=1}^{p^{-2}2^p}\pi_k \approx \frac{4A}{p}(\log(p)-\log\log(2p)-1)\hspace{15pt} A=\frac{p}{4(\log(p)-\log\log(2p)-1)}$
The renormalization is bizarre and I frankly don't trust it.  For one, it completely disregards the distinction between $2^p-1$, $(4^p-1)/3$, and $(16^p-1)/15$, since that boils down to the $4$ that gets swallowed by $A$.  For another, one gets $f(p) = O(p)$ which is much larger than observed data.
I think the conclusion of this, is that sieve methods are not applicable because the set of numbers we are considering do not constitute an appropriately dense portion of the integers.  It is not correct to give a sieve-like form $\approx 1/p_k$for the probability that $p_k$ divides our target number.  I'll keep thinking, but I'm posting my work here for people to at least learn a dead end.
