Integer-valued power sums

The function $$ f : z \in \mathbb{C} \longmapsto \sum_{i} \frac{a_i}{1-a_iz} $$ is meromorphic on $\mathbb{C}$ and has integral Taylor coefficients. It follows from a theorem of Borel that such a function must be in $\mathbb{Q}(z)$; see for example Richard Stanley's answer here. In particular $f$ has only finitely many poles ; this implies that only finitely many of the $a_i$'s are nonzero.

OK. I've finally got the time to correct my deleted answer. If anyone doesn't like complex analysis, this answer only uses real analysis.

Let $f(x)=\prod_{n=1}^\infty (1+a_nx)=\sum_{n=0}^\infty e_nx^n$, where $e_n$ is the elementary symmetric polynomial in $(a_n)$. Then as observed in the comments, $e_n\in\mathbb{Z}/n!$, so if infinitely many $a_n>0$, then $e_n\ge1/n!$, which implies $f(x)\ge e^x$, or $\log f(x)\ge x$ for $x>0$. On the other hand, we can pick $N$ such that $\sum_{n=N}^\infty a_n<1/2$. Then $\sum_{n=N}^\infty \log(1+a_nx)<x/2$. Also $\sum_{n=1}^{N-1} \log(1+a_nx)=o(x)$, so $\log f(x)\le(1/2+o(1))x$, contradiction.

This is not an answer --- but perhaps someone can fill the gap?

Clearly, we have $d>1$. Denote $$ n_r=\sum_{i=1}^\infty a_i^r\in\mathbb N. $$ Choose some $k$ such that $$ \sum_{i>k}a_i<\frac1{2d} $$ (all $k$ greater than some $k_0$ fit). Set $$ P(x)=(x-a_1)\dots(x-a_k)=x^k+p_{k-1}x^{k-1}+\dots+p_0. $$

Now, look at the rows $m_0,\dots,m_{3k}$ of the matrix $M_{k,s}=[n_{s+i+j}]_{0\leq i,j\leq 3k}$ (for some large $s$). We see that for every $i\geq k$ the elements of $$ m_i'=m_i+p_{k-1}m_{i-1}+\dots+p_0m_{i-k} $$ are expressed via the sums of large powers of $a_i$, $i>k$, with bounded coefficients. Thus all elements of $m_i'$ are less than $C_k/(2d)^s$ with some constant $C_k$ depending only on $k$. On the other hand, the entries of $m_0,\dots,m_{k-1}$ do not exceed $\alpha_k kd^{s+k}$ for some constant $\alpha_k$ also depending on $k$ only.

Replacing the $m_i$ with $m_i'$ does not chande $\det M_{s,k}$, so $$ \det M_{s,k}\leq (3k)!\cdot (\alpha_kkd^{s+k})^k\cdot \left(\frac{C_k}{(2d)^s}\right)^{2k} $$ which tends to $0$ as $s\to\infty$. On the other hand, $\det M_{s,k}$ is integer, so $\det M_{s,k}=0$ for all sufficiently large $s\geq s_k$.

Now I would like to say that this yields $(n_r)$ be a linear recurrent sequence from some moment. Ptifully, this is not true in general, but perhaps some argument may fill this?

Notice that if $(n_r)$ were linear recurrent, the rest is easy. Indeed, we may prove inductively (assuming $a_1\geq a_2\geq \dots$) that all nonzero $a_i$ are the roots of its characteristic polynomial, so they are finitely many.