  An improved measurement of the decay $B^0_s\to\mu^+\mu^-$ and searches for
the decays $B^0\to\mu^+\mu^-$ and $B^0_s\to\mu^+\mu^-\gamma$ are performed at
the LHCb experiment using data collected in proton-proton collisions at
$\sqrt{s} = 7,~8$ and $13$ TeV, corresponding to integrated luminosities of 1,
2 and 6 fb$^{-1}$, respectively. The $B^0_s\to\mu^+\mu^-$ branching fraction
and effective lifetime are measured to be
${\mathcal{B}}(B^0_s\to\mu^+\mu^-)=\left(3.09^{+0.46+0.15}_{-0.43-0.11}\right)\times
10^{-9}$ and $\tau(B^0_s\to\mu^+\mu^-)=(2.07\pm 0.29\pm 0.03)$ ps,
respectively, where the uncertainties include both statistical and systematic
contributions. No significant signal for $B^0\to\mu^+\mu^-$ and
$B^0_s\to\mu^+\mu^-\gamma$ decays is found and the upper limits
$\mathcal{B}(B^0\to\mu^+\mu^-)<2.6\times 10^{-10}$ and
$\mathcal{B}(B^0_s\to\mu^+\mu^-\gamma)<2.0\times 10^{-9}$ at 95% confidence
level are determined, where the latter is limited to the range $m_{\mu\mu} >
4.9$ GeV$/c^2$. Additionally, the ratio between the $B^0\to\mu^+\mu^-$ and
$B^0_s\to\mu^+\mu^-$ branching fractions is measured to be
$\mathcal{R}_{\mu^+\mu^-}<0.095$ at 95% confidence level. The results are in
agreement with the Standard Model predictions.
