A study of the decays B s0 → μ + μ − 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 s0 -B 0 mass difference, a single fit determines the signal yields for both decay modes. This results in a measurement of the branching fraction ℬ(Bs0→μ+μ−)=(3.2−1.0+1.1)×10−9 and an upper limit ℬ (B→ μ + μ − ) < 4.3 × 10 − 10 at 95% confidence level. The result is combined with the Run 1 ATLAS result, yielding ℬ(Bs0→μ+μ−)=(2.8−0.7+0.8)×10−9 and ℬ (B→ μ + μ − ) < 2.1 × 10 − 10 at 95% confidence level. The combined result is consistent with the Standard Model prediction within 2.4 standard deviations in the ℬ(B0→μ+μ−)−ℬ(Bs0→μ+μ−) plane.[Figure not available: see fulltext.]