Speed of convergence in Lebesgue's density theorem

The answer is negative. In fact, let us show that $a_n$ as defined in the question may converge to $0$ however slowly. Indeed, let $(\epsilon_j)$ is any sequence in $(0,1]$ converging to $0$ slowly and regularly enough, in the sense that \begin{equation} \epsilon_{j-1}-\epsilon_j\ge2^{2-j}\tag{1} \end{equation} eventually in $j$ -- that is, for some natural $j_0$ and all natural $j\ge j_0$. For instance, this condition will hold if $\epsilon_j$ is (eventually in $j$) of the form $2^{2-j}$ or $j^{-a}$ or $(\underbrace{\ln\dots\ln}_N j)^{-a}$ for some real $a$ and some natural $N$.

In what follows, $j$ is a natural number $\ge j_0$.
Let $$k_j:=\lfloor2^j\epsilon_j\rfloor$$ and $$r_j:=k_j/2^j.$$ Then \begin{equation} \epsilon_j-\tfrac1{2^j}<r_j\le\epsilon_j\tag{2} \end{equation} and hence, by the regularity condition $(1)$, $r_j<r_{j-1}$. Let $$\Delta_j:=(r_{j},r_{j-1}].$$ Note that the intervals $\Delta_j$ are pairwise disjoint. Let then $$s_{j,\ell}:=r_{j}+\ell/2^{j+1}$$ for $\ell=0,\dots,\ell_j:=(r_{j-1}-r_{j})2^{j+1}$, so that $\ell_j$ is even, and $$\delta_{j,\ell}:=(s_{j,\ell},s_{j,\ell+1}].$$ Note that $\Delta_j$ is the disjoint union of the $\delta_{j,\ell}$'s over $\ell=0,\dots,\ell_j-1$. Let, finally,
\begin{equation} B:=\bigcup_{j=j_0}^\infty\bigcup_{q=0}^{\ell_j/2-1}\delta_{j,2q}. \end{equation}

Take now any large enough natural $n$. Let \begin{equation} H_{n,i}:=\Big(\frac i{2^n},\frac{i+1}{2^n}\Big] \end{equation} for $i=0,\dots,2^n-1$.

From now on, suppose also that $j\ge n\ge j_0$.

Then $|\delta_{j,\ell}|=1/2^{j+1}<1/2^n=|H_{n,i}|$, where $|\cdot|$ denotes the Lebesgue measure. For any given $i=0,\dots,2^n-1$ such that $H_{n,i}\subseteq\Delta_j$, we have $$B\cap H_{n,i}=\bigcup_{q=q_1}^{q_2}\delta_{j,2q} \quad\text{and}\quad H_{n,i}=\bigcup_{\ell=2q_1}^{2q_2+1}\delta_{j,\ell}$$ for some integers $q_1=q_{1;n,i}$ and $q_2=q_{2;n,i}$ such that $0\le q_1\le q_2\le\ell_j/2-1$. Hence, $|B\cap H_{n,i}|=\frac12|H_{n,i}|=\frac12\frac1{2^n}$.
So, \begin{equation} \frac{|B\cap H_{n,i}|}{2^{-n}}\Big(1-\frac{|B\cap H_{n,i}|}{2^{-n}}\Big) =\frac14 \end{equation} -- for each $j\in\{n,n+1,\dots\}$ and each $i=0,\dots,2^n-1$ such that $H_{n,i}\subseteq\Delta_j$.
By $(2)$ and $(1)$, the length $|\Delta_j|$ of the dyadic interval $\Delta_j$ is $r_{j-1}-r_j\ge\epsilon_{j-1}-\epsilon_j-\tfrac1{2^{j-1}}\ge\tfrac12\,(\epsilon_{j-1}-\epsilon_j)$, and so, for a given natural $n\ge j_0$, the dyadic interval $\Delta_j$ contains at least $2^{n-1}(\epsilon_{j-1}-\epsilon_j)$ dyadic intervals $H_{n,i}$. So, for $a_n$ as defined in the question, we have
\begin{equation} a_n\ge\frac1{2^n}\sum_{j=n}^\infty\ \sum_{i\colon H_{n,i}\subseteq\Delta_j} \frac{|B\cap H_{n,i}|}{2^{-n}}\Big(1-\frac{|B\cap H_{n,i}|}{2^{-n}}\Big) \ge \frac1{2^n}\sum_{j=n}^\infty 2^{n-1}(\epsilon_{j-1}-\epsilon_j)\frac14 =\frac{\epsilon_{n-1}}8. \end{equation} Thus, $a_n\ge\frac{\epsilon_{n-1}}8$ for all $n\ge j_0$. For instance, letting $\epsilon_j=1/j$ for all natural $j$, we have
$\sum_n a_n=\infty$.

As Iosif Pinelis has just written, the answer is, indeed, negative. However, I would still like to make a couple of comments about this question which are too long for a comment.

First of all, this question is not about the Lebesgue density theorem. An important feature of this theorem is that it deals with neighborhoods of (almost) all points, and because of that it requires some information about the geometry of the space under consideration. However, the neighborhoods in the question are just the dyadic intervals which form a monotone sequence of partitions of the unit interval, so that the question belongs entirely to the measure category and is just about the properties of this sequence of partitions. In fact, the numbers $a_n$ are nothing else than $$ a_n = \| 1_B - \mathbf E(1_B|\mathfrak A_n) \|^2 \;, $$ where $\|\cdot\|$ is the $L^2$ norm on the unit interval, $\mathfrak A_n$ is the algebra of order $n$ dyadic intervals, and $f\mapsto \mathbf E(f|\mathfrak A_n)$ is the projection of $L^2$ onto the space of $\mathfrak A_n$-measurable functions given by the conditional expectation ($\equiv$ averaging over the order $n$ dyadic intervals). The convergence $$ \mathbf E(1_B|\mathfrak A_n) \to 1_B $$ in various senses (in particular, in $L^2$) is then just the theorem on convergence of conditional probabilities (more generally also known under the name of the martingale convergence theorem).

The example of Iosif Pinelis - if I understand it correctly - then boils down just to the following. First of all, there is no need to stick to the unit interval, and it is more convenient to describe it in (the equivalent) terms of the product (Bernoulli) measure space whose base is the two-point set $X=\{0,1\}$ endowed with the uniform measure. Then $\mathfrak A_n$ is just the algebra generated by the first $n$ coordinates. Now, let $A_n$ be a sequence of subsets of $X^n$, let $A'_n\subset X^{n+1}$ be the sets of strings obtained from $A_n$ by adding 1 in the $(n+1)$-th position, let $C'_n\in\mathfrak A_{n+1}$ be the associated cylinder sets, and let $$ B = \bigcup_n C'_n \;. $$ If the cylinder sets $C_n\in\mathfrak A_n$ associated with the sets $A_n$ are pairwise disjoint, then $|1_B-E(1_B|\mathfrak A_n)|$ vanishes on $C'_k$ for $k\le n-1$ and is 1/2 otherwise, whence $$ \sum \| 1_B - \mathbf E(1_B|\mathfrak A_n) \|^2 = \infty $$ iff $$ \sum n |C_n| = \infty \;. $$ These two conditions (disjointness of the sets $C_n$ and divergence of the above series) can be easily satisfied.