Asymptotics of an oscillatory integral with a linear oscillator I am interested in asymptotic results for
$$
S(p) = \int_0^1 \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} \sin(py) dy,
$$
i.e. a result that is valid as $p\rightarrow\infty$. The parameter $\varepsilon$ can be a complex number with $|\varepsilon|>1$ in the third quadrant of the complex plane. If you want a concrete number for plotting, take $\varepsilon=3$. I expect something like
$$
S(p) \approx \frac{1}{p^n}\sin(p+\alpha),
$$
for some $n$ and $\alpha$, both possibly functions of $\varepsilon$, but I don't know how to approach this. I know that this could be the leading order term of some kind of asymptotic series but don't know the details. Other similar integrals I've seen rely on the method of stationary phase, but those have some non-linear function of $y$ in the oscillatory part, and the minima of that function contribute to the integral. Can anyone enlighten me or point me to some references in the case of a simple sinusoidal oscillator?

Because the interval of integration is finite we can tackle this integral with the method of steepest descent.  We'll write $\DeclareMathOperator{rre}{Re} \DeclareMathOperator{iim}{Im}$
$$
\begin{align}
S(p) &= \int_0^1 \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} \sin(py) dy \\
&= \frac{1}{2i} \left( \int_0^1 \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} e^{ipy} dy - \int_0^1 \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} e^{-ipy} dy \right) \\
&= \frac{1}{2i} \Bigl( I_1(p) - I_2(p) \Bigr) \tag{1}
\end{align}
$$
and study $I_1$ and $I_2$ independently.  (If $\varepsilon^2$ were real we would only need to consider $\iim I_1(p)$, but let's study both to be safe.)
These integrals exist as long as $\varepsilon^2$ does not lie in the real interval $(-\infty,0)$.  Let $\pm y_\varepsilon$ be the roots of the equation
$$
(\varepsilon^2-1)y^2+1 = 0
$$
(assuming $\varepsilon \neq \pm 1$) and define the real number
$$
a = a(\varepsilon) = \begin{cases}
1 & \text{if } -1 \leq \varepsilon \leq 1, \\
\tfrac{1}{2} |\iim y_\varepsilon| & \text{otherwise.}
\end{cases}
$$
By definition we have $a > 0$.
Let's put $I_1$ into more standard terms by writing
$$
I_1(p) = \int_0^1 f(y) e^{pg(y)}\,dy.
$$
The current integration contour runs along a path of constant altitude on the surface described by
$$
\rre g(y) = \rre(iy) = -\iim y,
$$
and at the endpoints of the contour (namely $y=0$ and $y=1$) the paths of steepest descent are rays parallel to the positive imaginary axis.  We therefore deform the contour according to the following picture, where the original contour (the segment from $y=0$ to $y=1$) is shown in blue:

By our choice of $a$ we avoid the poles of the integrand and may thus conclude that
$$
I_1(p) = \int_0^1 = \int_{\gamma_1} + \int_{\gamma_2} - \int_{\gamma_3}. \tag{2}
$$
The contours $\gamma_1$ and $\gamma_3$ follow the paths of steepest descent, and the integral over the contour $\gamma_2$ is exponentially small:
$$
\left|\int_{\gamma_2} \frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1} e^{ipy} dy \right| \leq e^{-ap} \int_{\gamma_2} \left|\frac{y \sqrt{1-y^2}}{(\varepsilon^2-1)y^2+1}\right| |dy|. \tag{3}
$$
We can parameterize the contour $\gamma_1$ by $y = it$, $dy = idt$, to get
$$
\int_{\gamma_1} = -\int_0^a \frac{t \sqrt{1+t^2}}{1-(\varepsilon^2-1)t^2} e^{-pt}\,dt,
$$
which can be handled using Watson's lemma.  Indeed,
$$
\int_{\gamma_1} = -\frac{1}{p^2} + O\left(\frac{1}{p^4}\right). \tag{4}
$$
Similarly, $\gamma_3$ can be parameterized by $y = 1+it$, $dy = idt$, to get
$$
\int_{\gamma_3} = e^{i(p+\pi/2)} \int_0^a \frac{(1+it)\sqrt{1-(1+it)^2}}{(\varepsilon^2-1)(1+it)^2 + 1} e^{-pt}\,dt,
$$
which can again be handled using Watson's lemma.  We calculate
$$
\int_{\gamma_3} = e^{i(p+\pi/4)} \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} + O\left(\frac{1}{p^{5/2}}\right). \tag{5}
$$
Combining $(3)$, $(4)$, and $(5)$ in $(2)$ we find that
$$
I_1(p) = e^{i(p-3\pi/4)} \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} - \frac{1}{p^2} + O\left(\frac{1}{p^{5/2}}\right). \tag{6}
$$
An analogous method can be applied to $I_2$ using the contour in the lower half-plane joining the points $0, -ia, 1-ia, 1$.  Doing this, we obtain the asymptotic
$$
I_2(p) = e^{-i(p-3\pi/4)} \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} - \frac{1}{p^2} + O\left(\frac{1}{p^{5/2}}\right). \tag{7}
$$
Substituting $(6)$ and $(7)$ into $(1)$ we conclude that

$$
S(p) = \sin\left(p - \frac{3\pi}{4}\right) \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} + O\left(\frac{1}{p^{5/2}}\right). \tag{8}
$$

Here is a plot of the numerical evaluation of $S(p)$ in blue and the asymptotic in $(8)$ in red.  In this we've taken $\varepsilon = 3$.

The asymptotic $(8)$ is slightly out of phase with $S(p)$.  This is corrected in the next term of the asymptotic expansion, given by

$$
S(p) = \sin\left(p - \frac{3\pi}{4}\right) \frac{\sqrt{\pi/2}}{\varepsilon^2 p^{3/2}} + \cos\left(p - \frac{3\pi}{4}\right) \frac{3(8-3\varepsilon^2)\sqrt{\pi/2}}{8 \varepsilon^4 p^{5/2}} + O\left(\frac{1}{p^{7/2}}\right).
$$
  $\tag{9}$

