  We derive the self-consistent harmonic approximation for the $2D$ XY model
with non-local interactions. The resulting equation for the variational
couplings holds for any form of the spin-spin coupling as well as for any
dimension. Our analysis is then specialized to power-law couplings decaying
with the distance $r$ as $\propto 1/r^{2+\sigma}$ in order to investigate the
robustness, at finite $\sigma$, of the Berezinskii-Kosterlitz-Thouless (BKT)
transition, which occurs in the short-range limit $\sigma \to \infty$. We
propose an ansatz for the functional form of the variational couplings and show
that for any $\sigma>2$ the BKT mechanism occurs. The present investigation
provides an upper bound for the lower critical threshold $\sigma^\ast=2$, above
which the traditional BKT transition persists in spite of the LR couplings.
