Squarefree Parts of Mersenne Numbers
I can't prove that the natural density is $0$ but I wanted to mention some well known facts which might seem to hold promise for a proof.
The idea is that perhaps we can restrict the possible prime divisors of $b_n$ enough to force $a_n$ to grow too fast to have positive density (say faster than $\frac{2^n}{n^2}$.)
$7 \mid M_n$ exactly if $3 \mid n$ but $7^2 \mid M_n$ exactly if $3\cdot 7=21 \mid n$. We might record this information as $[7,3,3 \cdot 7].$ Then for the first few primes the information is $$[3,2, 2\cdot 3],[5, 4 , 4\cdot 5], [7, 3 , 3\cdot 7], [11, 10 , 10\cdot 11], [13, 12 , 12\cdot 13],$$$$ [17, 8 , 8\cdot 17], [19, 18 , 18\cdot 19], [23, 11 , 11\cdot 23], [29, 28 , 28\cdot 29], [31, 5 , 5\cdot 31]$$
So for these primes, a necessary (but not sufficient) condition for $p^2$ to divide $M_n$ is that $p$ divide $n$. Is this true for all primes $p \ge 3$? No there are least the two exceptions $1093$ and $3511$ whose data is $[1093,364,364],[3511,1755,1755]$, however those are the only exceptions at least as far as $1.2 \times 10^{17}.$ This means that the prime divisors of $b_n$ form a subset of the prime divisors of $n$ (at least for primes up to that bound with the two mentioned exceptions.) That would make it hard for $b_n$ to be large compared to $M_n.$
In slightly more detail:
Let $\ell=\mathrel{ord}_m(2)$ be the smallest positive integer with $M_{\ell}$ a multiple of $m.$ Then $\ell$ is a divisor of $\varphi(m)$, the number of integers $1 \le i \le m-1$ relatively prime to $m,$ then, $M_n$ is a multiple of $m$ exactly when $n$ is a multiple of $\mathrel{ord}_m(2)$.
The primes $p=1093$ and $p=3511$ are Wieferich primes in that $\mathrel{ord}_p(2)=\mathrel{ord}_{p^2}(2).$ It is known (according to Wikipedia) that there are no other Wieferich primes at least as far as $1.2 \times 10^{17}.$ Some people believe that the number of Wieferich primes up to $N$ grows like $\log (\log (N)).$ As far as I know, there is no proof which definitively rules out that there might be only those two nor any that rules out the possibility that all but finitely many primes are Wieferich primes .
Either $\mathrel{ord}_{p^{e+1}}(2)=p \cdot\mathrel{ord}_{p^{e}}(2)$ or $\mathrel{ord}_{p^{e+1}}(2)= \mathrel{ord}_{p^{e}}(2)$ Only in the cases above of $e=1$ and $p=1093$ or $3511$ is it known that $\mathrel{ord}_{p^{e+1}}(2)= \mathrel{ord}_{p^{e}}(2)$
Further details:
For $p=1093$ if $p \mid M_n$ also $p^2 \mid M_n$ , this happens when $364 \mid n.$ However it is still possible to have $a_n$ divisible by $1093$: If $n$ is a multiple of $364\cdot 1093$ but not of $364 \cdot 1093^2$ then $M_n$ is a multiple of $1093^3$ but not of $1093^4.$
For each integer $n \ge 1$ the cyclotomic polynomials $\Phi_n(x)$ is monic of degree $\varphi(n)$ and irreducible (over $\mathbb{Z}$). These polynomials are implicitly defined by $x^n-1=\prod_{d|n}\Phi_d(x).$ Hence
$$M_n=\prod_{d|n}\Phi_n(2).$$
This means that if $\ell=\mathrel{ord}_m(n)$ then not only does $m \mid M_{\ell}$, in fact $m \mid \Phi_{\ell}(2)$
$\gcd(M_a,M_b)=M_{\gcd(a,b)}$
So the only time that $M_q$ might be prime is when $q$ is prime. However most Mersenne numbers $M_q$ with prime index are not Mersenne primes. However in all known cases $M_q$ is square free when $q$ is prime (an exception would require a prime index $q$ and a divisor $p |M_q$ which is a Wieferich prime. Since neither of $364,1755$ is one less than a prime this would be a new, as yet unknown, Wieferich prime.
For $q$ qrime, every divisor of $M_q$ is of the form $2qk+1.$
For factorizations of $M_n$ up to $n=1199$ one can go to the 2-Minus tables of the Cunningham Project.
Here's a simple proof using user43383's idea and a recent result of Andrew Granville: http://arxiv.org/abs/1212.6306
According to Theorem 1 of that paper, $2^n-1$ always has a primitive prime factor $p$ that occurs to an odd exponent, except when $n=1$ or $n=6$ (where there are no primitive prime factors at all). Primitive here means that $2$ has order $n$ modulo $p$. In particular, $p\equiv 1\pmod{n}$. Clearly, $p \mid a_n$. So if $q_n$ denotes the smallest prime congruent to $1$ modulo $n$, then $a_n \geq q_n$ for every $n \neq 1, 6$. And in fact, by a direct check, $a_6 = 7 = q_6$, so $a_n \geq q_n$ for every $n > 1$.
Trivially, $q_n \geq n+1$ for every $n$. In fact, it is usually much larger. The simple fact that the primes have density zero implies that for every positive integer $K$, one has $q_n > Kn$ apart from a set of $n$ of density zero. (To see this, note that if this inequality fails, then one of $n+1$, $2n+1$, $\dots$, or $(K-1)n+1$ is prime, and each of these conditions puts $n$ in a set of density zero.)
Hence: $a_n \geq n+1$ for all $n > 1$, and for every fixed $K$, $a_n > Kn$ for all $n$ outside of a set of density zero. These two facts are enough to imply that $\{a_n\}$ itself has density zero.
(Proof of the last bit: Fix $K$. Given a large $x$, if $a_n \leq x$, then the inequality $a_n \geq n+1$ forces $n \leq x$. If $n \leq x$ and $a_n \leq x$, then either $n < x/K$ or $a_n \leq Kn$. The former holds for at most $x/K$ values of $n$, and the latter holds only for $o(x)$ values of $n$, as $x\to\infty$. Thus, the number of distinct $a_n$ with $a_n \leq x$ is at most $x(1/K+o(1))$; so the upper density of $\{a_n\}$ is at most $1/K$. But this holds for all $K$.)
I think I can prove the bound $a_n> \prod 2^{\alpha_i}p_i^{\frac{(\alpha_i)(\alpha_i+1)}{2}}$ if $n=\prod p_i^{\alpha_i}$. It is certainly very far of the real $a_n$, but it is more than enough to prove $\delta_A=0$.
First a lemma:
$\frac{2^{p^d}-1}{2^{p^{d-1}}-1}$ is never a square
Proof: If $p=2$ it is obviously true. Wlog $p$ is odd. We actually want to analyze the polynomial $f(x)=\frac{x^p-1}{x-1}$ for $x=2^{p^{d-1}}$. We will prove actually that this polynomial is never a square for "big" $x$. The idea is approximating $\sqrt{f(x)}$ by a polynomal with rational coefficients. Notice that $\sqrt{f(x)}<\sqrt{\frac{x^p}{x-1}}=x^{\frac{p-1}{2}}\sqrt{\frac{x}{x-1}}=x^{\frac{p-1}{2}}\sum\limits_{k=1}^{\infty} x^{-k}\frac{\binom{2k}{k}}{4^k}$ (by the generalized binomal theorem).
So $\sqrt{f(x)}-x^{\frac{p-1}{2}}\sum\limits_{k=1}^{\frac{p-1}{2}} x^{-k}\frac{\binom{2k}{k}}{4^k} < x^{\frac{p-1}{2}}\sum\limits_{k=\frac{p+1}{2}}^{\infty} x^{-k}\frac{\binom{2k}{k}}{4^k}<\frac{\binom{p+1}{\frac{p+1}{2}}}{2^{p+1}}\frac{1}{x-1}$ (substituting all denominators for the first one, which is greater).
Also, the left side can't be zero, because all the coefficients of $(x^{\frac{p-1}{2}}\sum\limits_{k=1}^{\frac{p-1}{2}} x^{-k}\frac{\binom{2k}{k}}{4^k})^2$ are at most $1$ (the first ones are exactly $1$, while the last ones are strictly less than $1$).
So if $\sqrt{f(x)}$ is an integer, the left side will be at least $\frac{1}{2^{p-1}}$ (because all denominators divide $2^{p-1}$). So, if $f(x)$ is a perfect square, $x-1<\frac{\binom{p+1}{\frac{p+1}{2}}}{2}<2^{p}-1$.
It indeed doesn't happen for $x=2^{p^{d-1}}$ for $d>1$. For $x=2$, $f(2)$ is clearly not a square because $f(2) \equiv 3$ mod $4$.
Now back to the bound.
Let $p^d$ a prime power factor of $n$. Notice that $gcd (\frac{2^{p^d}-1}{2^{p^{d-1}}-1},2^{p^{d-1}}-1)$ must divide $p$, (in general $gcd(\frac{x^p-1}{x-1},x-1)|p$, because $\frac{x^p-1}{x-1}=x^{p-1}+x^{p-2}+...+1 \equiv p \,(mod \,(x-1))$) therefore this gcd is $1$ (because clearly $p \nmid 2^{p^{d-1}}-1$).
Now we take a prime $q$ with odd exponent in $\frac{2^{p^d}-1}{2^{p^{d-1}}-1}$ (such prime exists because this expression is not a square). Since $q$ doesn't divide $2^{p^{d-1}}-1$ (by the earlier $gcd$ condition), $q$ satisfies $ord_q(2)=p^d \Rightarrow 2p^d|q-1 \Rightarrow q>2p^d$. Notice that this prime $q$ is a factor of $a_{p^d}$.
Repeating for $d-1$, $d-2$, ..., $1$ we get that $2^{p^d}-1$ has at least $d$ distinct prime factors that divide it with odd exponent, whose product is at least $2^d p^{\frac{d(d+1)}{2}}$. Since $2^{p_i^{\alpha_i}}-1|2^n-1$ and all $2^{p_i^{\alpha_i}}-1$ are pairwise coprime, we proved the bound previously stated. (notice that the factor $2^{\alpha_1}$ actually doesn't appear for $i=1$ ($p_1=2$), but it won't make a big difference).
Now, for the proof that $\delta_A=0$, we will use that the bound we proved gives $a_n>n 2^{\omega(n)-1}$. This is because $2^{\alpha_i}p_i^{\frac{(\alpha_i)(\alpha_i+1)}{2}}\ge 2 p_i^{\alpha_i}$ (for $i=1$, since we don't have the factor $p_i^{\alpha_i}$ we use only $p_1^{\frac{(\alpha_1)(\alpha_1+1)}{2}}\ge p_1^{\alpha_1}$, what explains the $-1$ after $\omega(n)$)
Now, set some fixed $t$:
\delta_A=$\lim\limits_{n \to \infty} \frac{\#\{ a_k \in A | a_k<n\}}{n} = \lim\limits_{n \to \infty} \frac{\#\{a_k \in A | a_k<n, \omega(k)<t\}}{n} + \frac{\#\{a_k \in A | a_k<n, \omega(k)\ge t\}}{n}$
By http://en.wikipedia.org/wiki/Hardy%E2%80%93Ramanujan_theorem, $\frac{\#\{a_k \in A | a_k<n, \omega(k)<t\}}{n}$ goes to zero, because by our bound, if $a_k<n$, then $k<n$, and one of the consequences of the theorem presented is that the set of numbers with $\omega(k)<t$ has natural density zero.
By our bound $\#\{a \in A | a<n, \omega(a)\ge t\}$ is at most $\frac{n}{2^{t-1}}$, because it has no elements $a_k$ with $k>\frac{n}{2^{t-1}}$ otherwise $a_k>k2^{\omega(k)-1}>\frac{n}{2^{t-1}}2^{t-1}=n $.
Therefore for each fixed $t$, $\delta_A\le\frac{1}{2^{t-1}}$, so we prove that indeed $\delta_A=0$.
(notice that the factor $2^{\alpha_i}$ in $2^{\alpha_i}p_i^{\frac{(\alpha_i)(\alpha_i+1)}{2}}$ is very important, otherwise we wouldn't be able to handle squarefree $n$)