I am looking for a simple proof of the following identity which holds by numerical tests:
$$ \int_{-\infty}^\infty x\left(\sqrt{1+\frac{a}{1-x^2-2i\gamma x}}-\sqrt{1+\frac{a}{1-x^2+2i\gamma x} }\right)dx={i\pi a}\operatorname{sgn}\gamma, $$ where $a$ and $\gamma$ are real numbers. The result suggests using the residue theorem, but the residues at singularities of the integrand seem to be zero...
Any hint or suggestion are appreciated.
The result depends on the choice of branch cut. Indeed, with the principal square root, the result is false if $a$ is negative with large modulus:
In this answer, we will prove a corrected version of OP's identity. To this end, define
$$ f(z) := 1+\frac{a}{1-z^2-2i\gamma z} \tag{1}. $$
This functions has poles at $z = -i\gamma \pm \sqrt{1-\gamma^2}$ and zeros at $z = -i\gamma \pm \sqrt{1+a-\gamma^2}$. Now if $a > -1$ and $\gamma > 0$, then all of these points lie in the lower half plane, and so, there exists a unique analytic function $g$ on the closed upper half plane $\overline{\mathbb{H}} = \{z \in \mathbb{C} : \operatorname{Im}(z) \geq 0\}$ such that
$$ \lim_{\overline{\mathbb{H}} \ni z \to \infty} g(z) = 1 \qquad \text{and} \qquad g(z)^2 = f(z). \tag{2} $$
Moreover, we can check that $g(x) = \sqrt{f(x)}$ holds for $x \in \mathbb{R}$, where $\sqrt{\cdot}$ is the principal square root. Now we claim
Note that this claim is enough to establish OP's identity for $a > -1$.
Proof of Proposition. Write $I$ for the left-hand side of $\text{(3)}$. Then
\begin{align*} I &= \lim_{R\to\infty} \int_{-R}^{R} x \left( \sqrt{f(x)} - \sqrt{f(-x)} \right) \, \mathrm{d}x \\ &= 2 \lim_{R\to\infty} \int_{-R}^{R} x \left( \sqrt{f(x)} - 1 \right) \, \mathrm{d}x \\ &= 2 \lim_{R\to\infty} \int_{-\Gamma_{R}} z \left( g(z) - 1 \right) \, \mathrm{d}z, \end{align*}
where $-\Gamma_{R}$ is the clockwise-oriented upper semicircular contour from $-R$ to $R$. Parametrizing $-\Gamma_{R}$ by $z = Re^{i\theta}$ with $\theta$ from $\pi$ to $0$, we get
\begin{align*} I &= -2i \lim_{R\to\infty} \int_{0}^{\pi} R^2e^{2i\theta} \left( g(Re^{i\theta}) - 1 \right) \, \mathrm{d}\theta. \end{align*}
But in view of $\text{(2)}$, it follows that
$$ R^2e^{2i\theta} \left( g(Re^{i\theta}) - 1 \right) = \frac{R^2e^{2i\theta} \left( f(Re^{i\theta}) - 1 \right)}{g(Re^{i\theta}) + 1} \xrightarrow[R\to\infty]{} -\frac{a}{2} $$
holds uniformly in $\theta \in [0, \pi]$. Therefore we can interchange the order of limit and integration to obtain
\begin{align*} I &= -2i \int_{0}^{\pi} \lim_{R\to\infty} R^2e^{2i\theta} \left( g(Re^{i\theta}) - 1 \right) \, \mathrm{d}\theta \\ &= -2i \int_{0}^{\pi} \left( -\frac{a}{2} \right) \, \mathrm{d}\theta \\ &= a\pi i. \end{align*}
Addendum. By the same technique, we have
The extra integral on the right-hand side comes from part of the branch cut $[0,i\beta]$ lying inside the upper half plane. And Mathematica seems to suggest that this integral reduces to elliptic integral.