Im trying to figure out a puzzling problem I've run into related to the series expansion for the Bessel function exponentiated, $(J_{v}(x))^{r}$ where $v \in \Bbb C$ and $r \in \Bbb R$. Using Euler's formula for such a thing, if $f(x) = \sum_{n=0}^{\infty} c_{n}x^{n}$ then $(f(x))^{r} = \sum_{n=0}^{\infty} d_{n}x^{n}$ where...
$$
d_{0} = c_{0}^{r} \quad , \quad d_{n} = \frac{1}{nc_{0}}\sum_{k=1}^{n} (k(r+1)-n)c_{k}d_{n-k}
$$
Then using this formula I found an expression for the Bessel function exponentiated to be...
$$
(J_{v}(x))^{r} = \frac{1}{(\Gamma(v+1))^{r}}\bigg(\frac{x}{2} \bigg)^{rv}\sum^{\infty}_{n=0} (-1)^{n}d_{n}\bigg(\frac{x}{2} \bigg)^{2n}
$$
$$
d_{0} = 1 \quad , \quad d_{n} = \frac{\Gamma(v+1)}{n}\sum_{k=1}^{n} \frac{(k(r+1)-n)}{k!\Gamma(k+v+1)}d_{n-k}
$$
My problem is that the terms of this series don't seem to be monotonically decreasing in computer implementations of the series. The behavior is such that the terms grow initially, which is to be expected, then the $d_{n}$ coefficients dominate the $x$ term then they bottom out at some value that depends on $x$, $r$, and $v$. Then each term continues to grow exponentially in $n$ again. The plots below demonstrate this.
This behavior doesn't seem to break having tested all the way up to $n=10,000$. I have implemented the code in three different languages (Julia, Python, and C++) and all have the same result and the formula doesn't seem to be wrong as it converges in a error bounds that is expected to the actual value if you take the terms before the minimum.
My question is what is wrong? Did I do the math wrong? Is there some condition that I'm missing in using Euler's formula above? Is there some nuance in computer implementations that Im missing? I've spent a week trying to figure out the problem and have been unsuccessful so any help would be much appreciated. Thus far I've tried to write a proof on the monotonicity of the terms but have fallen short and cant think what else to do. Thank you!