  A search for the fully reconstructed $B_s^0 \rightarrow \mu^+\mu^-\gamma$
decay is performed at the LHCb experiment using proton-proton collisions at
$\sqrt{s}=13$\,TeV corresponding to an integrated luminosity of
$5.4\,\mathrm{fb^{-1}}$. No significant signal is found and upper limits on the
branching fraction in intervals of the dimuon mass are set
  \begin{align}
  {\cal B}(B_s^0 \rightarrow \mu^+\mu^-\gamma) <
4.2\times10^{-8},~&m(\mu\mu)\in[2m_\mu,~1.70]\,\mathrm{GeV/c^2} ,\nonumber
  {\cal B}(B_s^0 \rightarrow \mu^+\mu^-\gamma) <
7.7\times10^{-8},~&m(\mu\mu)\in[1.70,~2.88]\,\mathrm{GeV/c^2},\nonumber
  {\cal B}(B_s^0 \rightarrow \mu^+\mu^-\gamma) <
4.2\times10^{-8},~&m(\mu\mu)\in[3.92 ,~m_{B_s^0}]\,\mathrm{GeV/c^2},\nonumber
\end{align} at 95\% confidence level. Additionally, upper limits are set on the
branching fraction in the $[2m_\mu,~1.70]\,\mathrm{GeV/c^2}$ dimuon mass region
excluding the contribution from the intermediate $\phi(1020)$ meson, and in the
region combining all dimuon-mass intervals.
