Recently I was learning some properties of asymptotic expansions. One property is: if $$f(z)\sim a_0+\frac{a_1}{z}+\frac{a_2}{z^2}+\cdots$$ as $z\to \infty$, then $$\int_z^\infty f(\zeta)-a_0-\frac{a_1}{\zeta} \,\mathrm{d}\zeta\sim \sum_{n=2}^\infty \frac{a_n}{(n-1)z^{n-1}}$$ holds independently of the path from $z$ to $\infty$ as long as all paths lie in a wedge-shaped region $\alpha<\arg(z)<\beta$.
The definition of asymptotic expansion is $$f(z)\sim a_0+\frac{a_1}{z}+\frac{a_2}{z^2}+\cdots$$ as $z\to \infty$ (The series is not necessarily convergent) iff for any $N$, $$\frac{f(z)-a_0-\cdots-\frac{a_N}{z^N}}{\frac{1}{z^N}}\to 0$$ as $z\to\infty$.
My attempt
$$|\int_z^\infty \left(f(\zeta)-a_0-\frac{a_1}{\zeta} \right)\mathrm{d}\zeta-\frac{a_2}{z}-\cdots-\frac{a_{N+1}}{N z^N}|= |\int_z^\infty \left(f(\zeta)-a_0-\frac{a_1}{\zeta}-\frac{a_2}{\zeta^2}-\cdots-\frac{a_{N+1}}{\zeta^{N+1}} \right)\mathrm{d}\zeta|$$ Fix $\epsilon>0$, it suffices to show the above expression is less than $\epsilon |z|^{-N}$ when $|z|$ large enough.
By definition of asymptotic expansions, we can find $R$ such that $$|f(z)-a_0-\frac{a_1}{z}-\frac{a_2}{z}-\cdots-\frac{a_{N+1}}{z^{N+1}}|<\epsilon |z|^{-N-1}$$ when $|z|>R$.
If I estimate the integrand by the bound, I get something like $$\int_z^\infty \epsilon |\zeta|^{-N-1} |\mathrm{d}\zeta|$$ If the path is along the real axis, the proof is finished. However, I have no idea how to deal with the case involved a complex path. Additionally, I'm not sure if $f(z)$ is analytic or if $\int_z^\infty f(\zeta)-a_0-\frac{a_1}{\zeta} \mathrm{d}\zeta$ is well defined(finite?, independent of path?) either.