  A study of the decays $B^0_s \to \mu^+\mu^-$ and $B^0 \to \mu^+\mu^-$ has
been performed using 26.3 fb$^{-1}$ of 13 TeV LHC proton-proton collision data
collected with the ATLAS detector in 2015 and 2016. Since the detector
resolution in $\mu^+\mu^-$ invariant mass is comparable to the
$B_{s}^{0}$-$B^0$ mass difference, a single fit determines the signal yields
for both decay modes. This results in a measurement of the branching fraction
${\cal B}(B^0_s \to \mu^+\mu^-) = \left( 3.2^{+1.1}_{-1.0} \right) \times
10^{-9}$ and an upper limit ${\cal B}(B^0 \to \mu^+\mu^-) < 4.3 \times
10^{-10}$ at 95% confidence level. The result is combined with the Run 1 ATLAS
result, yielding ${\cal B}(B^0_s \to \mu^+\mu^-) = \left( 2.8^{+0.8}_{-0.7}
\right) \times 10^{-9}$ and ${\cal B}(B^0 \to \mu^+\mu^-) < 2.1 \times
10^{-10}$ at 95% confidence level. The combined result is consistent with the
Standard Model prediction within 2.4 standard deviations in the ${\cal B}(B^0
\to \mu^+\mu^-)$-${\cal B}(B^0_s \to \mu^+\mu^-)$ plane.
