The sum of reciprocals of even index Lucas numbers has a nice closed-form in terms of theta functions, $$\begin{aligned}S_e &= \sum_{n=1}^\infty \frac1{L_{2n}}\\ &= \sum_{n=1}^\infty \frac1{\phi^{2n}+\phi^{-2n}}\\ &= \tfrac14\Big(\vartheta_3^2(\phi^{-2})-1\Big)\\ &= 0.56617\dots \end{aligned}$$
with golden ratio $\phi$ and Jacobi theta function $\vartheta_3(q)$.
However, it seems this is just a special case. Given the root $x_1>0$ and $x_2<0$ of the quadratic, $$x^2+bx-1=0$$
Is it true that, $$\sum_{n=1}^\infty \frac1{x_1^{2kn}+x_2^{2kn}} = \tfrac14\Big(\vartheta_3^2(x_1^{-2k})-1\Big) $$ for any integer $k>0$, and where the Lucas numbers was just $b=-1$?
Since $x_1^{2k}x_2^{2k} = (-1)^{2k} = 1$, we may choose $q = x_i^{2k}$ such that $|q| \leq 1$. (For Lucas numbers we have $q = x_2^{2k} = x_1^{-2k}$. But I guess this choice would depend on $b$.) Since $|q| = 1$ case makes the sum diverge, let us focus on the case where $|q| < 1$. Then it boils down to showing
$$ \sum_{n=1}^{\infty} \frac{q^n}{1+q^{2n}} = \frac{1}{4}(\vartheta_3(q)^2 - 1 ). $$
It is well-known that $ r_2(n) = \#\{ (a,b)\in\mathbb{Z}^2 : a^2+b^2=n\}$ satisfies
$$ \forall n \geq 1 \ : \quad r_2(n) = 4 \sum_{\substack{d \mid n \\ d \text{ odd}}} (-1)^{\frac{d-1}{2}}. $$
So it follows that for $|q| < 1$,
\begin{align*} \vartheta_3(q)^2 &= \sum_{(a,b)\in\mathbb{Z}^2} q^{m^2+n^2}=\sum_{n=0}^{\infty} r_2(n) q^n \\ &= 1 + 4\sum_{j=1}^{\infty}\sum_{k=0}^{\infty} (-1)^k q^{j(2k+1)} \\ &= 1 + 4\sum_{j=1}^{\infty} \frac{q^j}{1+q^{2j}}. \end{align*}