I claim, for $artan(x) = \frac{1}{2}ln(\frac{1+x}{1-x})$, i.e., the inverse hyperbolic tangent function:
$$
y=1+b_1artan(x)+b_2artan(x)^2=1+a_1x+a_2x^2+a_3x^3 +...a_jx^j...
$$

for $a_j > 0$. In other words, the given quadratic is expressible exactly as a power series. Assuming this is so, how can I relate $(b_1, b_2)$ to $a_j$, so that knowing the one allows me to compute the other? If this is not solvable for the given information, is it solvable in the special case of the power series being a geometric series, so that all $a_j$ are equal? I'm aware there is a relationship between Taylor series and power series that may be exploitable to express a polynomial as a power series, but not sure how or whether it relates here.

The question arises in the theory of nonuniform random permutations. Let $X$ be a permutation of the first $n$ nonzero positive integers, let $Y$ be a random permutation of $X$, and let their expected correlation be $-1 <\rho_{XY}<1$. Given certain distributional assumptions, it is known that, for infinite $n$, the expected number of fixed points is the sum of the geometric series in $\rho_{XY}$ with coefficients $a = 1$. For finite $n$, it is known the expected number of fixed points is the sum of a power series in $\rho_{XY}$ with all coefficients nonnegative.

To my knowledge, the expression for the power series coefficients in the finite case is an open problem, although clearly each would have to be a function of $n$ that goes to 1 as $n$ goes to infinity. If my conjecture holds, then for $x = \rho_{XY}$, I can approximate the polynomial coefficients via Monte Carlo simulation. But I still need to convert them to approximate power series coefficients, which is where you come in!