I can find the principal part of the Laurent expansion $$ \frac{1}{(J_1 (z-r_k))^2}=\frac{a_{-2}^k}{(z-r_k)^2}+\frac{a_{-1}^k}{z-r_k}+\sum_{j=0}^\infty a_j^k(z-r_k)^j $$ where $r_k$ is one of the simple zeros of $J_1(z)$ without trouble using the power series for $J_1(z)$. With these I can define the "analytic remainder" $$ h(z)=\frac{1}{(J_1 (z-r_k))^2}-\frac{a_{-2}^k}{(z-r_k)^2}-\frac{a_{-1}^k}{z-r_k}=\sum_{j=0}^\infty a_j^k(z-r_k)^j $$ Estimating $h(z)$ as $|z|\rightarrow\infty$ is easy because then $h(z)\rightarrow\frac{1}{(J_1 (z-r_k))^2}$, as the principal part becomes small. (Of course $\frac{1}{(J_1 (z-r_k))^2}$ has infinitely many arbitrarily large and real poles $r_k$, but the point here is that estimation in this case is not a problem.)
However, I am not sure how to estimate $h(z)$ in the neighborhood of $r_k$. I could in principle use power series to calculate the other coefficients $a_j^k$ but this seems a very inefficient and inelegant way of doing so. What would be a more effective way?