How do I solve the following definite integral (preferably by an asymptotic method)? 
$$ \int_{1}^{\infty} \frac{\sin^2 (\mu \sqrt{x^2 -1})}{(x+1)^{\frac{9}{2}} (x-1)^{\frac{3}{2}}} \,dx $$
Note: $\mu$ here is an extremely small constant.

I have tried:

*

*Estimating the integral by Taylor expansion of $\sin^2(\mu \sqrt(x^2 - 1)$ but the it diverges after few terms.

*I have also tried to split the integral into two regions, solving the two integrals differently and then adding them in the hope that the arbitrary parameter $\alpha$ will vanish (asymptotic matching/splitting). While I have been able to solve the two integrals but $\alpha$ does not vanish: $$\int_{1}^{\alpha} f(x)\,dx + \int_{\alpha}^{\infty} f(x)\,dx$$
Approximations in for the two regions:
$1$ to $\alpha$ region: $x \tilde = 1$ and hence $(x+1) \tilde= 2$ 
$\alpha$ to $\infty$ region: $(x^2 -1) \tilde= x^2$
 A: Here is a log-log plot of
$$\delta I=I_{\text{appr}}-\int_{1}^{\infty} \frac{\sin^2 (\mu \sqrt{x^2 -1})}{(x+1)^{\frac{9}{2}} (x-1)^{\frac{3}{2}}} \,dx,$$
as a function of $\mu$, with
$$I_{\text{appr}}=\frac{2 \mu^2}{15}-\frac{\mu^4}{9}.$$
You write "$\mu$ is extremely small". Is this error acceptable?

This is a plot of the absolute error; the relative error is $\lesssim 10^{-6}$ for $\mu\lesssim 10^{-2}$.
A: It is more or less straightforward to write down an asymptotic expansion around $\mu\to0$,
$$
I(\mu)\sim \frac{2 \mu ^2}{15}-\frac{\mu ^4}{9}+\frac{\pi  \mu ^5}{15}+\frac{1}{450} \mu ^6 (60 \log (\mu )+60 \gamma -67)+\cdots
$$
where $\gamma$ is Euler-Mascheroni. Given this result it should be rather clear why the naive Taylor expansion diverges after the first two terms: the next terms have odd powers of $\mu$ and even non-analytic factors of $\log(\mu)$, while your function $f$ only has even powers of $\mu$.
Here is a log-plot of the error:

(Blue means one-term approximation, Red is two-term, Green is three-term, and Orange is four-term.)
A simple Mathematica code to reproduce the formula above (and easily generalized to yield more terms):
Sin[µ Sqrt[x^2 - 1]]^2/((x + 1)^(9/2) (x - 1)^(3/2))
Series[%, {µ, 0, 10}]
small = Integrate[Normal[%], {x, 1, Λ}, Assumptions -> Λ > 1];
Sin[µ x]^2 Series[1/((x + 1)^(9/2) (x - 1)^(3/2)), {x, ∞, 8}] // Normal
large = Integrate[%, {x, Λ, ∞}, Assumptions -> Λ > 1 && 0 < µ < 1];
Assuming[Λ > 1 && 0 < µ < 1, Series[small + large, {µ, 0, 10}] // Normal] // Short
Assuming[Λ > 1 && 0 < µ < 1, Series[%, {Λ, ∞, 1}] // Normal]

A: I really like the double-series approach by @AccidentalFourierTransform, however I got remaining $\Lambda$s in the terms of order $O(\mu^8)$ onwards. Thinking about this approach, I located the problem in the neglection of higher order terms in $\sin^2(\mu x)$ in the series expansion in $x$ around infinity (line 4 in the MMA code, thanks for providing!).
The problem can be eliminated by a change of variables in the original integral from $x \mapsto y \equiv \sqrt{x^2-1}$, such that
$$
I = \int_0^\infty \mathrm dy \frac{y \sin^2(\mu y)}
{x\,(x-1)^{3/2}\,(x+1)^{9/2}}, \quad \text{with} \quad x \equiv \sqrt{y^2+1}.\tag{1}
$$
Now, the method can be applied: we split the integral at $\Lambda>0$, expand the first integrand in $\mu$ around 0 and the second integrand in $x$ around $\infty$. The sum of the integrated series expansions is again expanded around $\mu=0$ to get
\begin{align}
I&=\frac{2\mu^2}{15} - \frac{\mu^4}{9} + \frac{\pi\mu^5}{15}
- \frac{[67 - 60(\gamma+\log\mu)]\,\mu^6}{450}
- \frac{8\pi\mu^7}{315}\\
&\quad + \frac{[169 - 56 (\gamma+\log\mu)]\,\mu^8}{7056}
+ \frac{[913 - 360 (\gamma+\log\mu)]\,\mu^{10}}{2916000} \\
&\quad - \frac{[23797 - 9240 (\gamma+\log\mu)]\,\mu^{12}}{3841992000}
+ O\left(\mu ^{14}\right).\tag{2}
\end{align}
Note that now $\Lambda$ cancels in all terms (as suggested by the OP in 2.) and that the highest odd power seems to be $\mu^7$.
