A study of the decays B s 0 → μ + μ − and B 0 → μ + μ − 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 μ + μ − 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 ℬ B s 0 → μ + μ − = 3.2 − 1.0 + 1.1 × 10 − 9 $$ \mathrm{\mathcal{B}}\left({B}_s^0\to {\mu}^{+}{\mu}^{-}\right)=\left({3.2}_{-1.0}^{+1.1}\right)\times {10}^{-9} $$ and an upper limit ℬ B 0 → μ + μ − < 4.3 × 10 − 10 $$ \mathrm{\mathcal{B}}\left({B}^0\to {\mu}^{+}{\mu}^{-}\right)<4.3\times {10}^{-10} $$ at 95% confidence level. The result is combined with the Run 1 ATLAS result, yielding ℬ B s 0 → μ + μ − = 2.8 − 0.7 + 0.8 × 10 − 9 $$ \mathrm{\mathcal{B}}\left({B}_s^0\to {\mu}^{+}{\mu}^{-}\right)=\left({2.8}_{-0.7}^{+0.8}\right)\times {10}^{-9} $$ and ℬ B 0 → μ + μ − < 2.1 × 10 − 10 $$ \mathrm{\mathcal{B}}\left({B}^0\to {\mu}^{+}{\mu}^{-}\right)<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 ℬ B 0 → μ + μ − − ℬ B s 0 → μ + μ − $$ \mathrm{\mathcal{B}}\left({B}^0\to {\mu}^{+}{\mu}^{-}\right)-\mathrm{\mathcal{B}}\left({B}_s^0\to {\mu}^{+}{\mu}^{-}\right) $$ plane.