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 Wierferich primes in that $\mathrel{ord}_p(2)=\mathrel{ord}_{p^2}(2).$ It is known (according to Wikipedia) that there are no other Wierferich primes at least as far as $1.2 \times 10^{17}.$ Some people believe that the number of Wierferich 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 Wierferich 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 Wierferich prime. Since neither of $364,1755$ is one less than a prime this would be a new, as yet unknown, Wierferich 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.

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.

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}$.) I think what follows is worth mentioning. However , if I understand things correctly (see below) it is not known that $A$ is infinite. That certainly argues against being able to prove that $a_n$ grows.

Let $\ell=\mathrel{ord}_m(2)$ be the smallest positive integer with $M_{\ell}$ a multiple of $m.$ Then $\ell$ is a divisor of $\phi(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)$.

I mentioned the extremely unlikely, but not, to my knowledge, disproved possibility that only finitely many primes fail to be Wierferich primes. If I reason correctly, that would make $A$ have density $0$ because, rather than growing, it is finite.

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}$.)

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)$.

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 Wierferich prime. Since neither of $364,1755$ is one less than a prime this would be a new, as yet unknown, Wierferich 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.