How can I simplify this sum any further?

You might be able to use the fact that $$\sum_{k=0}^\infty b_{ak}=\sum_{k=0}^\infty \left(\frac{1}{a}\sum_{j=0}^{a-1} \exp\left(2\pi ijk/a\right)\right)b_k.$$ For example, when $a=1$, taking $b_k = \frac{x^k}{k!}\sum_{\ell \ge 0} \binom{k}{\ell}$ yields $$\sum_{k=0}^\infty b_{k}=\sum_{k=0}^\infty \frac{x^k}{k!}\sum_{\ell \ge 0} \binom{k}{\ell}=\sum_{k=0}^\infty \frac{x^k}{k!}2^k=\exp(2x),$$ as you already obtained. For $a=2$, first note that \begin{align} \sum_{\ell \ge 0}\binom{k}{2 \ell} &= \sum_{\ell\ge 0} \left(\frac{1}{2}\sum_{j=0}^1 \exp\left(2\pi ij\ell/2\right)\right)\binom{k}{\ell}\\ &= \sum_{\ell\ge 0} \frac{1+(-1)^\ell}{2}\binom{k}{\ell}\\ &= \frac{1}{2}\sum_{\ell\ge 0} \binom{k}{\ell}+ \frac{1}{2}\sum_{\ell\ge 0} (-1)^\ell\binom{k}{\ell}\\ &= \frac{2^k+0^k}{2}. \end{align} Now taking $b_k = \frac{x^k}{k!}\sum_{\ell \ge 0}\binom{k}{2 \ell}$ yields \begin{align}\sum_{k=0}^\infty b_{2k}&=\sum_{k=0}^\infty \left( \frac{1}{2}\sum_{j=0}^1 \exp\left(\pi ijk\right)\right)\frac{x^k}{k!}\sum_{\ell \ge 0}\binom{k}{2 \ell}\\ &=\frac{1}{2}\sum_{k=0}^\infty \left(1+(-1)^k\right)\frac{x^k}{k!}\left(2^{k-1}+\frac{1}{2}[k=0]\right)\\ &=\frac{1}{4}\sum_{k=0}^\infty \frac{x^k}{k!}2^k+\frac{1}{4}\sum_{k=0}^\infty (-1)^k\frac{x^k}{k!}2^k+\frac{1}{2}\\ &=\frac{\exp(2x)+\exp(-2x)}{4} +\frac{1}{2}\\ &=\cosh^2(x), \end{align} again matching your result. For $a=3$, first note that \begin{align} \sum_{\ell \ge 0}\binom{k}{3 \ell} &= \sum_{\ell\ge 0} \left(\frac{1}{3}\sum_{j=0}^2 \exp\left(2\pi ij\ell/3\right)\right)\binom{k}{\ell}\\ &= \sum_{\ell\ge 0} \frac{1+\exp(2\pi i\ell/3)+\exp(4\pi i\ell/3)}{3}\binom{k}{\ell}\\ &= \frac{1}{3}\sum_{\ell\ge 0} \binom{k}{\ell}+ \frac{1}{3}\sum_{\ell\ge 0} \exp(2\pi i/3)^\ell\binom{k}{\ell}+ \frac{1}{3}\sum_{\ell\ge 0} \exp(4\pi i/3)^\ell\binom{k}{\ell}\\ &= \frac{2^k+(1+\exp(2\pi i/3))^k+(1+\exp(4\pi i/3))^k}{3}\\ &= \frac{2^k+\exp(\pi i/3)^k+\exp(-\pi i/3)^k}{3}. \end{align} Now taking $b_k = \frac{x^k}{k!}\sum_{\ell \ge 0}\binom{k}{3 \ell}$ yields \begin{align}\sum_{k=0}^\infty b_{3k}&=\sum_{k=0}^\infty \left( \frac{1}{3}\sum_{j=0}^2 \exp\left(2\pi ijk/3\right)\right)\frac{x^k}{k!}\sum_{\ell \ge 0}\binom{k}{3 \ell}\\ &=\frac{1}{3}\sum_{k=0}^\infty \left(1+\exp(2\pi ik/3)+\exp(4\pi ik/3)\right)\frac{x^k}{k!}\frac{2^k+\exp(\pi i/3)^k+\exp(-\pi i/3)^k}{3}\\ &=\frac{1}{9}\sum_{k=0}^\infty (1+\exp(2\pi i/3)^k+\exp(4\pi i/3)^k)(2^k+\exp(\pi i/3)^k+\exp(-\pi i/3)^k)\frac{x^k}{k!}. \end{align} Now expand the product of trinomials to obtain 9 sums that reduce to $\exp(cx)$ for various constants $c$.

Alternatively, note that: $$\sum_{k=0}^\infty \frac{x^{ak}}{(ak)!}\sum_{\ell \ge 0}\binom{ak}{a\ell} = \left(\sum_{k=0}^\infty \frac{x^{ak}}{(ak)!}\right)^2,$$ so you might as well just compute \begin{align} \sum_{k=0}^\infty \frac{x^{ak}}{(ak)!} &= \sum_{k=0}^\infty \left( \frac{1}{a}\sum_{j=0}^{a-1} \exp\left(2\pi ijk/a\right)\right)\frac{x^k}{k!} \\ &= \frac{1}{a}\sum_{j=0}^{a-1} \sum_{k=0}^\infty \frac{(\exp(2\pi ij/a)x)^k}{k!} \\ &= \frac{1}{a}\sum_{j=0}^{a-1} \exp(\exp(2\pi ij/a)x), \end{align} and then square the result.


The explicit formula is as follows: $$ S_a=\frac{1}{a^2}\left(\sum_{z^a=2^a}+2\sum_{p_a(z)=0}\right)e^{az} $$ where the polynomials $p_a$ are given by A244608. For example, \begin{align} p_9(z)&= 1 - 13604 z^9 - 13359 z^{18} + 247 z^{27} + z^{36}\\ p_{10}(z)&=3125-383750 z^{10}-73749 z^{20}+502 z^{30}+z^{40} \end{align} whose roots in the complex plane look like follows:

enter image description here

enter image description here

The first few solutions are \begin{align} S_1&=e^{2 x}\\ S_2&=\frac{e^{-2 x}}{4}+\frac{e^{2 x}}{4}\\ S_3&=\frac{2 e^{-x}}{9}+\frac{e^{2 x}}{9}+\frac{2}{9} e^{\frac{x}{2}-\frac{1}{2} i \sqrt{3} x}+\frac{2}{9} e^{\frac{x}{2}+\frac{1}{2} i \sqrt{3} x}+\frac{1}{9} e^{-x-i \sqrt{3} x}+\frac{1}{9} e^{-x+i \sqrt{3} x}\\ S_4&=\frac{e^{-2 x}}{16}+\frac{1}{8} e^{(-1-i) x}+\frac{1}{8} e^{(-1+i) x}+\frac{1}{16} e^{-2 i x}+\frac{1}{16} e^{2 i x}+\frac{1}{8} e^{(1-i) x}+\frac{1}{8} e^{(1+i) x}+\frac{e^{2 x}}{16} \end{align} as given by \begin{align} p_1(z)&=0\\ p_2(z)&=0\\ p_3(z)&=1+z^3\\ p_4(z)&=4+z^4\\ p_5(z)&=-1+11z^5+z^{10}\\ p_6(z)&=-27+26z^6+z^{12} \end{align} etc. Quoting the OEIS entry, the coefficients are found as follows:

Let $\omega$ be a primitive $j$-th root of unity. Let $L(k)=\sum_{p=0}^{j-1} c(p)\omega^{kp}$ with $c(0)=2$ and $c(i)=C(j,i)$ if $i>0$. Then $p(j,X)=(X-L(1))(X-L(2))\dots(X-L([(n-1)/2]))$.


Actually the true generalisation for your $a=1$, $a=2$ cases is this expression for your sum: $$\frac{1}{a^2}\sum_{s=0}^{a-1}\sum_{r=0}^{a-1}e^{x\omega^r(1+\omega^s)}.$$ where $\omega$ is an $a$-th root of unity ($\omega^a=1$, $a \neq 1$).

When $a=1$ we have $\omega=1$ and your sum is $e^{2x}$. For $a=2$ we have $\omega=-1$ and your sum is $\frac{1}{4}(e^{2x}+e^{-2x}+e^0+e^0)=\frac{1}{4}(e^{2x}+e^{-2x}+2)$. Both agreeing with your observations.

Proof: Using the fact that $\sum_{l=0}^{a-1}\omega^{lr}=0$ unless $a|l$ we can show: $$\sum_{l=0}^k\binom{ak}{al}=\frac{1}{a}\sum_{s=0}^{a-1}(1+\omega^s)^{ak}$$ and$$\sum_{k=0}^\infty\frac{x^{ak}}{(ak)!}=\frac{1}{a}\sum_{r=0}^{a-1}e^{\omega^r x}.$$

Substituting the above expressions into your sum we obtain:

$$\sum_{k=0}^\infty\frac{x^{ak}}{(ak)!}\biggl(\sum_{l=0}^k\binom{ak}{al}\biggr)=\sum_{k=0}^\infty\frac{x^{ak}}{(ak)!}\frac{1}{a}\sum_{s=0}^{a-1}(1+\omega^s)^{ak}=\frac{1}{a}\sum_{s=0}^{a-1}\sum_{k=0}^\infty\frac{(x(1+\omega^s))^{ak}}{(ak)!}$$$$=\frac{1}{a}\sum_{s=0}^{a-1}(\frac{1}{a}\sum_{r=0}^{a-1}e^{x\omega^r(1+\omega^s)})=\frac{1}{a^2}\sum_{s=0}^{a-1}\sum_{r=0}^{a-1}e^{x\omega^r(1+\omega^s)}.$$