Let $f(x)=\dfrac{a^2}{x^2-a^2}$. I want to evaluate the Mellin transform of this function. Let me denote it ${\cal M}[f]$. What I want is to evaluate $${\cal M}[f](s)=\int_0^\infty dx x^{s-1}\dfrac{a^2}{x^2-a^2}.$$
Looking at the form of the integral and considering previous experience with Fourier transforms, I would consider evaluating this using contours in the complex plane. Still, here the integral runs over $(0,\infty)$ so I'm unsure how to proceed.
Suppose that $a\in \mathbb{R}$. Then the integrand has poles at $x=\pm a$. Having that in mind one possibility that I have considered was to take the contours (where $R > a$)
- $\Gamma_1$: parameterized by $t\mapsto t$ with $t\in [0,a-\epsilon]$;
- $\Gamma_2$: parameterized as $\theta\mapsto a+\epsilon e^{i\theta}$ where $\theta\in [-\pi,0]$;
- $\Gamma_3$: parameterized by $t\mapsto t$ where $t\in [a+\epsilon,R]$;
- $\Gamma_4$: parameterized by $\theta\mapsto Re^{i\theta}$ where $\theta\in [0,\frac{\pi}{2}]$;
- $\Gamma_5$: parameterized by $t\mapsto i(R-t)$ where $t\in [0,R]$;
Now let $\Gamma$ be obtained by joining all of these contours. There is no pole in the region enclosed by $\Gamma$ and so the integral of $g(z)=z^{s-1}\dfrac{a^2}{z^2-a^2}$ over $\Gamma$ is zero. As a result we have that when we send $R\to \infty$ and $\epsilon\to 0$ the integrals over $\Gamma_1$ and $\Gamma_3$ should sum up to ${\cal M}[f](s)$ and therefore we would have something of the form $${\cal M}[f](s)=-\lim_{R\to \infty}\lim_{\epsilon \to 0}\left(\int_{\Gamma_2}g(z)dz+\int_{\Gamma_4}g(z)dz+\int_{\Gamma_5}g(z)dz\right)$$
Still I'm unsure. This seems complicated at first. Is this really the way to evaluate ${\cal M}[f]$? Is there a better way of evaluating this Mellin transform?
I don't know if the contour integral approach will work. However, I want to suggest another way you can obtain the transform.
Assume $a\in\mathbb{R}\setminus\{0\}$ (Same as you assumed, but without 0 as @Maxim pointed out). Make the change of variables: $$ \begin{aligned} x&=|a|t^{1/2}\\ dx &= \frac{|a|}{2}t^{-1/2}dt \end{aligned} $$ In this sense, if $t\geq 0$ then $x\geq 0$ also. So that: $$ \begin{aligned} \int \frac{a^2x^{s-1}}{x^2-a^2}dx &= \int \frac{|a|^2x^{s-1}}{x^2-|a|^2}dx=\int \frac{x^{s-1}}{(x/|a|)^2-1}dx \\&= \frac{|a|^{s}}{2}\int t^{\frac{s-1}{2}}(t-1)^{-1}t^{-1/2}dt = -\frac{|a|^{s}}{2}\int t^{\frac{s}{2}-1}(1-t)^{-1}dt \end{aligned} $$ Now, recall the definition of the incomplete Beta function: $$ B(z;\alpha,\beta) = \int_0^zt^{\alpha-1}(1-t)^{\beta-1}dt $$ Thus, if we consider first the interval $[0,L]$ instead of the whole $[0,\infty)$ we get: $$ \int_0^L \frac{a^2x^{s-1}}{x^2-a^2}dx = -\frac{|a|^{s}}{2}\int_0^{(L/|a|)^2}t^{\frac{s}{2}-1}(t-1)^{-1}dt = -\frac{|a|^{s}}{2} B\left((L/|a|)^2;\frac{s}{2},0\right) $$ Now we just need to compute: $$ \lim_{L\to\infty} -\frac{|a|^{s}}{2} B\left((L/|a|)^2;\frac{s}{2},0\right) $$ To do so, use this identity from Wikipedia here (which is a simple change of variables): $$ B(z;\alpha,\beta) = (-1)^\alpha B\left(\frac{z}{z-1};\alpha,1-\alpha-\beta\right) $$ Applying it: $$ \begin{aligned} \lim_{L\to\infty} -\frac{|a|^{s}}{2} B\left((L/|a|)^2;\frac{s}{2},0\right) &= \lim_{L\to\infty}-\frac{|a|^{s}(-1)^{\frac{s}{2}}}{2} B\left(\frac{(L/|a|)^2}{(L/|a|)^2-1};\frac{s}{2},1-\frac{s}{2}\right)\\ &= -\frac{|a|^{s}(-1)^{\frac{s}{2}}}{2}B\left(1;\frac{s}{2},1-\frac{s}{2}\right) \end{aligned} $$ And, fortunately, $B(1;\alpha,\beta)=B(\alpha,\beta)$ with $B(\bullet,\bullet)$ the typical beta function: $$ B(1;\alpha,\beta)=B(\alpha,\beta):=\int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}dt = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} $$ where $\Gamma(\bullet)$ is the Gamma function. Thus, $$ -\frac{|a|^{s}(-1)^{\frac{s}{2}}}{2}B\left(1;\frac{s}{2},1-\frac{s}{2}\right) = -\frac{|a|^{s}(-1)^{\frac{s}{2}}}{2}\Gamma\left(\frac{s}{2}\right)\Gamma\left(1-\frac{s}{2}\right) $$ since $\Gamma\left(\frac{s}{2} + 1-\frac{s}{2}\right) = \Gamma(1)=1$. Now, we can further simply by applying the identity in section 18.4.5 here: $$ \Gamma(z)\Gamma(1-z) = \pi\text{csc}(z) $$ So that: $$ -\frac{|a|^{s}(-1)^{\frac{s}{2}}}{2}\Gamma\left(\frac{s}{2}\right)\Gamma\left(1-\frac{s}{2}\right) = -(-1)^{\frac{s}{2}}\frac{\pi |a|^{s}}{2}\text{csc}\left(\frac{\pi s}{2}\right) $$ And we can conclude: $$ \mathcal{M}\left[\frac{a^2}{x^2-a^2}\right](s) = -(-1)^{\frac{s}{2}}\frac{\pi |a|^{s}}{2}\text{csc}\left(\frac{\pi s}{2}\right) $$ This coincides with the result from Wolfram alpha here.
Hope this helps!