Given $a \ge 0$ and $b \gt 1$, and large $n$, how to approximate $x$; $1 =\frac{1}{n^ax}+\frac{1}{n}\sum_{k=1}^\infty\frac{1}{1+k^bx}$ Fix real numbers $a \ge 0$ and $b \gt 1$. For any integer $n\ge 1$, let $x_n$ be the unique positive solution of the equation
$$
\tag{1}
1 = \frac{1}{n^a x} + \frac{1}{n}\sum_{k=1}^\infty \frac{1}{1 + k^b x}.
$$
Of course, $x_n$ depends on $a$ and $b$. It is easy to see (check (1) below) that $x_n \ge n^{-a}$. Less trivial is the fact that
$$
\frac{x_n}{n^{-a}} - 1 = O(n^{-(1-a)}), \text{ if }a \in (0,1).
\tag{0}.
$$
Indeed, from (1), observe that
$$
x - n^{-a}= \frac{x}{n}\sum_{k=1}^\infty \frac{1}{1 + k^b x} =  \frac{1}{n }\sum_{k=1}^\infty  k^{-b} \frac{k^b x}{1 + k^b x} \le \frac{1}{n }\sum_{k=1}^\infty  k^{-b} =   O(1/n).
$$

Question. For large $n$, what is good and asymptotically exact formula for $x_n$ ? That is, construct $\epsilon_n$ such that $x_n/n^{-a} - 1 = O(\epsilon_n)$, with $\epsilon_n$ decays as fast as possible.

Note. Fixed-point equations like (1) appear in random matrix theory.
A possibly sub-optimal solution

Claim. If $b \in (1,\infty)$ and $a \in [0,b)$, then for large $n$, it holds that
$$
\frac{x_n}{n^{-a}} - 1=O(n^{- (1- a/b)}) = o(1).
$$

Proof. Indeed, since the second term in the RHS of (1)  is nonnegative, we deduce that
$$
x_n \ge n^{-a}.
\tag{2}
$$
Now, observe that because the function $t \mapsto 1 / (1 + xt^b)$ is decreasing on $[0,\infty)$, one has
$$
\sum_{k=1}^\infty \dfrac{1}{1 + k^b x}  \le \int_0^\infty \dfrac{1}{1 + z^b x}\mathrm{d}z  = cx^{-1/b},
$$
where $c:=\dfrac{\pi/b}{\sin(\pi/b)} \in (0,\infty)$, and the last step is thanks to this ME post https://math.stackexchange.com/a/247903/168758. Combining with (1) then gives
$$
1 = \frac{1}{n^a x_n} + O(\frac{1}{nx_n^{1/b}}).
\tag{3}
$$
Combining (2) and (3) gives $1 \le n^{-a}/x_n + O(n^{-(1-a/b)})$, and so
$$
\frac{x_n}{n^{-a}} \le \frac{1}{1-O(n^{-(1-a/b)})} 
 = 1+O(n^{-(1-a/b)}),
\tag{4}
$$
Combining with (1) gives the claimed result. $\quad\Box$

Points of improvement
I don't think I've made the best use of (3). Maybe it is possible to obtain a better $\epsilon_n$ via a more clever use of (3).
 A: I think that the idea of using the continuous case is a way to go for possible approximations. Multiplying by $x$ to remove the asymptote, we need to find the zero of function
$$f(x)= x-\frac{\pi   \csc \left(\frac{\pi }{b}\right)}{b\, n}x^{1-\frac{1}{b}}-n^{-a}$$ Let
$$x=\frac y{n^a}\qquad \qquad \alpha=\frac{\pi}b \csc \left(\frac{\pi }{b}\right)\,n^{\frac{a}{b}-1}\qquad \qquad \beta=1-\frac 1b$$
$$f(y)=y-\alpha \,y^\beta-1$$ The first derivative cancels at
$$y_*=(\alpha  \beta )^{\frac{1}{1-\beta }}$$ and then a first estimate using series
$$y_0=y_*+\sqrt{-2 \frac{f(y_*)}{f''(y_*)}}$$
and the first iterate of Newton method
$$y_1=y_0-\frac{f(y_0)}{f'(y_0)}$$ is totally explicit.
Trying with  $(a=\pi,b=e,n=11)$
$$y_*=1.49611\quad \quad  y_0= 5.39684 \quad \quad y_1=7.86537 \quad \implies \quad x_1=0.00420812$$ while the solution is $y=7.62709$ that is to say $x=0.00408063$.
Trying with the summation, the solution is $x=0.00369922$
Edit
Probably easier would be to let $$y=z^{\frac{1}{\beta }}\qquad \qquad \gamma=\frac{b}{b-1} >1$$
$$f(z)=z^\gamma-\alpha z-1$$
