Exact sum rule for solutions of the transcendental equation $x \tan x=\xi$

Question

Consider the transcendental equation \begin{equation} x_j \tan x_j=\xi \end{equation} with $\xi$ real and positive and with $j=0,1,2,\dots$. The roots lie approximately at \begin{equation} x_j \approx \begin{cases} &\sqrt{\frac{\xi}{3+\xi}}, \quad j=0\,\newline & \pi j + \frac{\xi}{\pi j}, \quad j=1,2,\dots \end{cases} \end{equation} see Eq. (13) of this article. For $\xi\leq 1$, the error of this approximation is less then $1\%$. The error increases as $\xi$ becomes larger.

I believe$^{1,2}$ that the roots satisfy the exact sum rule \begin{equation} \sum_{j=0}^{\infty}\frac{\sin^2 x_j}{\xi + \sin^2 x_j} =\frac{1}{2}. \end{equation} My question is whether there is a direct route (as opposed to footnote 1) to prove this sum rule. Are there other examples of solutions to transcendental equations satisfying exact sum rules?

$^1$ The above sum rule follows from the high frequency ($\omega\to\infty$) limit of the impedance of the so-called transmission line model, see this article. That article shows numerically in Fig. 10 that (I dropped some unimportant factors) \begin{equation} Z(i\omega)=\left( \sum_{j\ge0}\frac{4 x_j \sin^2 x_j}{2x_j +\sin 2 x_j} \frac{i \omega }{x_j^2 +i \omega }\right)^{-1} \end{equation} is identical to \begin{equation} Z(i\omega)=\frac{1}{\xi}+\sqrt{\frac{1}{i \omega} } \coth \sqrt{ i \omega } \end{equation} For the case $\xi\to\infty$, the equivalence of the two $Z(i\omega)$ expressions can be shown analytically, see Eqs. (61)-(64) of the mentioned article. I would be interested to see a proof for the equivalence of two $Z(i\omega)$ expressions for arbitrary $\xi$.

$^2$ The script below supports this conjecture.

import numpy as np
from scipy import optimize

def transcendental_eq(x, xi):
    return x*np.sin(x)-xi*np.cos(x)

def x_sol(xi):
    return [optimize.newton(transcendental_eq,i, args=(xi,)) for i in x0]

ml=1000 #sum range 
x0= np.concatenate(([np.pi/2], [np.pi*i+np.pi/4 for i in range(1,ml)]))#guess for solution

def sumrule(xi):
    x = x_sol(xi)
    s=0
    for i in range(0,ml):
        s+=np.sin(x[i])**2/(xi+np.sin(x[i])**2)
    return s    

print(sumrule(0.1))
print(sumrule(1))
print(sumrule(10))

Output:

0.4999898628139413

0.4998986281525962

0.49898628584402044