  We determine, by means of lattice QCD calculations, the local form factors
describing the $B_{s}\to \mu^{+}\mu^{-}\gamma$ decay. For this analysis we make
use of the gauge configurations produced by the ETM Collaboration with
$N_{f}=2+1+1$ flavour of Wilson-Clover twisted-mass fermions at maximal twist.
To obtain the $B_{s}$ meson form-factors, we perform simulations for several
heavy-strange meson masses $m_{H_{s}}$ in the range $m_{H_{s}} \in [ m_{D_{s}},
2 m_{D_{s}} ]$, and extrapolate to the physical $B_{s}$ meson point
$m_{B_{s}}\simeq 5.367~{\rm GeV}$ making use of the HQET scaling laws. We cover
the region of large di-muon invariant masses $\sqrt{q^{2}} > 4.16\,{\rm GeV}$,
and use our results to determine the branching fraction for $B_{s}\to
\mu^{+}\mu^{-}\gamma$, which has been recently measured by LHCb in the region
$\sqrt{q^{2}} > 4.9\,{\rm GeV}$. The largest contribution to the uncertainty in
the partial branching fractions at values of $\sqrt{q^{2}} < 4.8\,{\rm GeV}$ is
now due to resonance and other long-distance effects, including those from
"charming penguins", which we estimate by summing over the contributions from
the $J_P=1^-$ charmonium resonances.
