Experimenting with a CAS suggests an induction. In order to handle the induction, we need to consider the forms of the numbers involved. $\frac{4^m-1}{3} = 1 + 2^2 + 2^4 + \cdots + 2^{2m-2}$ alternates $1$ and $0$ bits. The map $2n + 1 \to n - 2^{f(n)}$ drops the rightmost bit (which is $1$) and clears the next least significant bit. The map $2n \to n$ drops the rightmost bit (which is $0$). And the map $2n \to 2n - 2^{f(n)}$ moves the least significant bit right one place. So the numbers which occur in the evaluation tree of $\frac{4^m-1}{3}$ are of the form $2^i + 2^j + 2^{j+2} + \cdots + 2^{j+2k}$ where $i \le j - 2$ and $k \ge 0$; and (when we get down to one bit) the form $2^i$.
Starting with the simplest case, when $i > 0$, $a(2^i) = pa(2^{i-1}) + qa(2^{i-1})$ so by induction $a(2^i) = (p+q)^i$.
For the more general case, let $A(i,j,k) = a(2^i + 2^j + 2^{j+2} + \cdots + 2^{j+2k})$ subject to the aforementioned constraints on $i,j,k$.
For $i > 0$, we have an even argument and use the second case of the recursion:
\begin{eqnarray*}
A(i,j,k) &=& a(2^i + 2^j + 2^{j+2} + \cdots + 2^{j+2k}) \\
&=& pa(2^{i-1} + 2^{j-1} + 2^{j+1} + \cdots + 2^{j+2k-1}) + qa(2^{i-1} + 2^j + 2^{j+2} + \cdots + 2^{j+2k}) \\
&=& pA(i-1, j-1, k) + qA(i-1, j, k)
\end{eqnarray*}
Then it's a trivial proof by induction that $$A(i,j,k) = \sum_{u=0}^i \binom{i}{u} p^u q^{i-u} A(0, j-u, k)$$
For $i = 0$, we have an odd argument and use the third case of the recursion:
\begin{eqnarray*}
A(0,j,k) &=& a(1 + 2^j + 2^{j+2} + \cdots + 2^{j+2k}) \\
&=& a(2^{j+1} + \cdots + 2^{j+1+2(k-1)}) \\
&=& \begin{cases}
a(0) & \textrm{if } k=0 \\
a(2^{j+1}) & \textrm{if } k=1 \\
A(j+1, j+3, k-2) & \textrm{otherwise}
\end{cases} \\
&=& \begin{cases}
1 & \textrm{if } k=0 \\
(p+q)^{j+1} & \textrm{if } k=1 \\
\sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} A(0, j+3-u, k-2) & \textrm{otherwise}
\end{cases}
\end{eqnarray*}
Further CAS experimentation suggests that the theorem we need to prove is $$A(0, j, 2v-1) = A(0, j, 2v) = \left(p\frac{q^v-1}{q-1} + q^v\right)^{j+1} A(0, 2, 2v-2)$$
The first of those equalities is easy: $$A(0,j,2) = \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} A(0, j+3-u, 0) = \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} = (p+q)^{j+1}$$ so $A(0,j,1) = A(0,j,2)$ and since the only occurrence of $k$ in the third case is in the parameter $k-2$ the rest follows by induction on $v$.
The second equality is the interesting one. Again, by induction on $v$:
\begin{eqnarray*}
A(0,j,2v) &=& \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} A(0, j+3-u, 2v-2) \\
&=& \sum_{u=0}^{j+1} \binom{j+1}{u} p^u q^{j+1-u} \left(p\frac{q^{v-1}-1}{q-1} + q^{v-1}\right)^{j+4-u} A(0, 2, 2v-4) \\
&=& \left[ \sum_{u=0}^{j+1} \binom{j+1}{u} p^u \left(pq\frac{q^{v-1}-1}{q-1} + q^v\right)^{j+1-u} \right] \color{Blue}{\left(p\frac{q^{v-1}-1}{q-1} + q^{v-1}\right)^3 A(0, 2, 2v-4)} \\
&=& \left( p + pq\frac{q^{v-1}-1}{q-1} + q^v \right)^{j+1} \color{Blue}{A(0, 2, 2v-2)} \\
&=& \left(p\frac{q^v-1}{q-1} + q^v\right)^{j+1} A(0, 2, 2v-2)
\end{eqnarray*}
as desired.
The answer to the original question now drops out: $a\left(\tfrac{4^m-1}{3}\right) = A(0, 2, m-2)$, but $A(0, 2, k)$ has the form of a cube times $A(0, 2, k-2)$ so by induction it's always a cube.
The conjectured formula can be proved by induction on $\mathrm{wt}(n)$.
For $\mathrm{wt}(n)=0$, we have $n=0$ and the conjectured formula trivially holds.
Now, for a positive integer $\ell$, suppose that the conjectured formula holds for all $\mathrm{wt}(n)<\ell$. Let us prove it for $\mathrm{wt}(n)=\ell$. So, let $n=2^{t_1}(1+2^{t_2+1}(1+\dots(1+2^{t_\ell+1}))\dots)$ with $t_j\geq 0$ and $n'=2^{t_2}(1+\dots(1+2^{t_\ell+1}))\dots)$, so that
\begin{split}
a(n) &= a(2^{t_1}(2n'+1)) = \sum_{m=0}^{t_1} \binom{t_1+1}{m} a(2^kn') \\
&=\sum_{m=0}^{t_1} \binom{t_1+1}{m} \sum_{j=0}^{2^{\ell-1}-1} (-1)^{\ell-1-\mathrm{wt}(j)} (1+\mathrm{wt}(j))^{m+t_{2}+1} \prod_{k=2}^{\ell} (1+\mathrm{wt}(\lfloor \tfrac{j}{2^{k-1}}\rfloor))^{t_{k+1}+1} \\
&=\sum_{j=0}^{2^{\ell-1}-1} (-1)^{\ell-1-\mathrm{wt}(j)} \left[(2+\mathrm{wt}(j))^{t_{1}+1} - (1+\mathrm{wt}(j))^{t_{1}+1}\right] \prod_{k=1}^{\ell} (1+\mathrm{wt}(\lfloor \tfrac{j}{2^{k-1}}\rfloor))^{t_{k+1}+1} \\
&=\sum_{j=0}^{2^{\ell-1}-1} \left[(-1)^{\ell-\mathrm{wt}(2j+1)}\prod_{k=0}^{\ell} (1+\mathrm{wt}(\lfloor \tfrac{2j+1}{2^k}\rfloor))^{t_{k+1}+1} + (-1)^{\ell-\mathrm{wt}(2j)} \prod_{k=0}^{\ell} (1+\mathrm{wt}(\lfloor \tfrac{2j}{2^k}\rfloor))^{t_{k+1}+1}\right] \\
&=\sum_{j'=0}^{2^\ell-1} (-1)^{\ell-\mathrm{wt}(j')} \prod_{k=0}^{\ell} (1+\mathrm{wt}(\lfloor \tfrac{j'}{2^k}\rfloor))^{t_{k+1}+1}, \\
\end{split}
where $j'$ corresponds to values $2j$ and $2j+1$ (in other words, $j=\lfloor j'/2\rfloor$). QED
Best Answer
The definition of $a_1$ given in OEIS is based on a bijection between integer partitions and natural numbers. A partition $\lambda_1\geq\lambda_2\geq\dots\geq\lambda_m>0$ with exactly $m$ parts corresponds to the number $$2^{\lambda_1+m-2}+2^{\lambda_2+m-1}+\dots+2^{\lambda_m-1}.$$ The definition of $a_1$ can then be written as \begin{equation}\label{a1d}(1)\qquad a_1\left(2^{\lambda_1+m-2}+2^{\lambda_2+m-1}+\dots+2^{\lambda_m-1}\right)=\lambda_1\dotsm\lambda_m.\end{equation}
When $n$ is a natural number, I will write $\|n\|$ for the number of ones in the binary expansion of $n$, and $n\preceq m$ if all binary digits in $n$ are smaller than or equal to those in $m$. It is well-known that $\binom nk\,\operatorname{mod}\,2=1$ if and only if $k\preceq n$. It is easy to see from this that $$a_m(n)=\sum_{k\preceq n}(m-1)^{\|n-k\|}a_1(k).$$ Roughly speaking, we go from $k$ to $n$ by changing zeroes to ones, and for each of the $\|n-k\|$ zeroes we can choose to change it in any one of $m-1$ binomial transforms.
The final ingredient is the identity $$\left\{\array{k+l\\l}\right\}=\sum_{l\geq \lambda_1\geq\lambda_2\geq\dots\geq\lambda_k> 0}\lambda_1\lambda_2\dotsm\lambda_k.$$ Probably this is well-known. I verified it using induction; just let me know if you need more explanation. By (1), this can be written \begin{equation}\label{sa}(2)\qquad\left\{\array{k+l\\l}\right\}=\sum_{0\leq j<2^{l+k-1},\,\|j\|=k}a_1(j).\end{equation}
We now have all ingredients we need. We have $$s_m(n)=\sum_{0\leq k<2^n}a_m(k)=\sum_{0\leq k<2^n}\sum_{l\preceq k}(m-1)^{\|k-l\|}a_1(l).$$ For fixed $l$ and $j=\|k-l\|$, we obtain $k$ by choosing $j$ from $n-\|l\|$ zeroes in $l$ and changing them to ones. Thus, we can write $$s_m(n)=\sum_{0\leq l<2^n}a_1(l)\sum_{j}\binom{n-\|l\|}{j}(m-1)^{j} =\sum_{0\leq l<2^n}a_1(l)m^{n-\|l\|}.$$ Writing this as a sum over $k=n-\|l\|$ and using (2) gives indeed $$s_m(n)=\sum_{k=0}^n m^k\left\{\array{n+1\\k+1}\right\}. $$