Is it true that the set $$S:=\{n\in \mathbb N\ |\ 2^n-1 \ \mbox{ is square free}\}$$ has positive density? What can we say when we replace $2^n-1$ with $\frac{a^n-1}{a-1}$?



Density of all n such that 2^n-1 is square free
User43383's answer linking the convergence of the sum of the reciprocals of the order of $2 \pmod{p^2}$ to the existence of a positive density is very nice. Let me add a little to this by saying what the density should be heuristically (approximately $0.754$), and by obtaining upper bounds for this density unconditionally (approximately $0.761$). For an odd prime $p$ let $\omega(p)$ denote the order of $2 \pmod{p^2}$. If $k$ is a square-free odd integer 
define $\omega(k)$ to be the l.c.m. of all the $\omega(p)$'s with $p|k$; thus $\omega(k)$ is the order of $2 \pmod{k^2}$. 


Note that $k^2$ divides $2^n-1$ if and only if $\omega(k)$ divides $n$. Now note that $\omega(3)=6$, $\omega(5)=20$, $\omega(7)=21$, $\omega(11)=110$, $\omega(13)=156$ and so on. Thus if $n$ is divisible by $6$, or $20$ or $21$, or $110$, or $156$ etc then $2^n-1$ is not squarefree. The density of the set of multiples of $\omega(p)$ for $3\le p\le z$ say may be calculated by using inclusion-exclusion. Put $P(z) = \prod\_{3\le p\le z} p$; then this density is 
$$ 
\sum\_{d|P(z)} \frac{\mu(d)}{\omega(d)}.
$$  

It is not immediate that this density is at least 
$$ 
\prod\_{3\le p\le z} \Big(1-\frac{1}{\omega(p)}\Big)
$$ 
as indicated in User43383's answer, but this is true by an inequality of Heilbronn and Rohrbach (see Theorem 0.9 in Hall's book Sets of Multiples; Cambridge tract 118). In any case, computing this density for any $z$ would give an upper bound for the density of squarefree Mersenne numbers. Thus taking $z=31$ one gets an upper bound of $0.761645...$ for the density. 


So far our argument is easy to make rigorous. Now for the heuristic part. The density computed above is a monotone decreasing function of $z$ clearly. Therefore it converges to some number, which will be strictly positive if $\sum 1/\omega(p)$ converges (thanks to the Heilbronn-Rohrbach inequality and the argument of User43383). One would expect that this limiting value is the right density. 


Alternatively we can argue as follows: Using that $\sum\_{k^2|n} \mu(k) = 1$ if $n$ is square-free and $0$ otherwise, we may 
write 
$$
\sum\_{n\le N} \mu(2^n-1)^2 = \sum\_{k \text{ odd}} \mu(k) \sum\_{n\le N, k^2|2^n-1} 1.
$$ 
The inner sum over $n$ may be expected to be roughly $N/\omega(k)$. Omitting error terms, which are significant since the sum over $k$ is over a large range, we may expect an asymptotic of 
$$ 
N \sum\_{k \text{ odd}} \frac{\mu(k)}{\omega(k)}.
$$ 
In other words, the density should be 
$$ 
\sum\_{k \text{ odd} } \frac{\mu(k)}{\omega(k)}, 
$$ 
which presumably converges, and presumably to the same answer obtained by considering $P(z)$ as before. 


**Numerics:** There are only two known Wieferich primes $1093$ and $3511$ up to $4 \times 10^{12}$. I calculated the conjectured density above up to $10000$ and the probability seems to be around $0.754$ (the first two digits seem reliable and the third fluctuates a bit). I also found $7550$ square-free Mersenne numbers with $n$ up to $10000$ (to be precise the numbers I found here are highly likely to be squarefree -- they are not divisible by the squares of the known Wieferich primes, and also not by the squares of all non Wieferich primes -- but there is a possibility that some of them could be divisible by the square of some ginormous Wieferich prime); the fit seems very good. 


**Remarks:** I also note that the convergence of the sum $\sum\_p 1/\omega(p)$ has come up earlier in work of Granville and Soundararajan on an Erdos problem of writing odd numbers as a squarefree plus a power of $2$. See <http://www.dms.umontreal.ca/~andrew/PDF/wieferich.pdf> for a copy of that paper, especially page 9 there which shows under reasonable conjectures on Wieferich primes that $\sum\_p 1/\omega(p)< 0.9091$. 


Let me also refer to a somewhat related MO problem [Probability of coprime polynomials](https://mathoverflow.net/questions/145757/probability-of-coprime-polynomials/146050#146050)
for which I discuss a heuristic similar to the one above. That problem may be viewed as a function field analog of asking what is the probability that $2^n-1$ and $3^m-1$ are coprime if $m$ and $n$ are chosen randomly. In problems of this type the expected densities do not appear to be products over primes as one may first think, since the divisibility conditions imposed by different primes are highly dependent.


