You do not state any conditions on the functions on which your operator is defined (it involves an infinite series which must converge). So my answer will also be purely formal.
Setting $g(t)=f(e^t)$, we obtain the equation
$$\sum_n a_n g(t+c_n)=\lambda g(t),$$
where $c_n=\log n$. This is a linear equation and such equations are solved with Fourier transform. Or Laplace transform whichever is more appropriate after you choose your space of functions. See, for example,
On equation f(z+1)-f(z)=f'(z)
You will get plenty of eigenfunctions. Indeed, taking FT, gives
$$(\sum_n a_n e^{ic_nz}-\lambda)\hat{g}(z)=0.$$
The expression in parentheses is an entire function which has infinitely many zeros. For every such zero $z_k$, we have a solution $\hat{g}(z)=\delta(z-z_k)$. The inverse Fourier transform of delta is an exponential. So the general solution is an exponential sum.
Once we know that, we can forget about Fourier, and just plug $g(t)=e^{pt}$ to the equation and determine $p$, as we do with linear ODE with constant coefficients in our undergraduate classes.