Density of all n such that 2^n-1 is square free 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}$?
 A: To my best  knowledge, we don't know even if the set $S$ has an infinity of elements.
In other words we don't know if there exist infinitely many squarefree numbers of the form $2^n-1$.
So,i think this is related to an open problem.
A: I give here an heuristic argument that this density should be positive. And we shouldn't expect any non heuristic argument for now, since as Lucia said, the infinitude of $S$ implies the existence of infinitely many non Wieferich primes.
Notice that $2^n-1$ is squarefree iff there is no prime $p$ such that $ord_{p^2}(2)|n$. Heuristically, the set of such numbers has a density of at least:
$ \prod\limits_{prime \,p} (1-\frac{1}{ord_{p^2}(2)})$ $\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(*)$
So it is enough to prove this infinite product is positive. Taking $log$ and using $log(1-x) \sim -x$ for small $x$, it is enough to prove that:
$\sum\limits_{prime \, p} \frac{1}{ord_{p^2}(2)}$ converges.
Now, notice that if $p$ is Wieferich prime, $ord_{p^2}(2)=ord_{p}(2)$, and if $p$ is a non Wieferich prime, $ord_{p^2}(2)=p\, ord_{p}(2)$. Notice also that $ord_p(2)>\log_2(p)$
Here we separate the sum between Wieferich primes and non Wieferich primes. For the first we get:
$\sum\limits_{Wieferich\, p}  \frac{1}{ord_{p^2}(2)}<
\sum\limits_{Wieferich\, p}  \frac{1}{log_2(p)}$
But by the classical heuristic that Wieferich primes should have density $\frac{log(log(n))}{n}$, the $n$-th Wieferich prime should be at least about $e^{e^n}$ $(**)$ , so $log_2(p)$ should grow exponentially, and the above sum should be geometric and converge. And for the second:
$\sum\limits_{non\, Wieferich\, p} \frac{1}{ord_{p^2}(2)}<
\sum\limits_{non\, Wieferich\, p} \frac{1}{p\,log_2(p)} \ll
\sum\limits_{n=2}^{\infty} \frac{1}{n\,log^2(n)}$
That converges by comparison with $\int \frac{1}{x\,log^2(x)}$
Note: Once we notice that $\sum\limits_{prime \, p} \frac{1}{ord_{p^2}(2)}$ converges, the step $(*)$ can be done non heuristically, so the heuristic part of this argument is only $(**)$.
A: 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
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.
