\noindent A search for the decays $B^0_{s}\rightarrow \mu^+ \mu^ \mu^+ \mu^$ and $B^0 \rightarrow \mu^+ \mu^ \mu^+ \mu^$ is performed using data, corresponding to an integrated luminosity of 1.0 $\{ fb}^{1}$ , collected with the LHCb detector in 2011. The number of candidates observed is consistent with the expected background and, assuming phasespace models of the decays, limits on the branching fractions are set: \{${\cal B }(B^0_{s}\rightarrow \mu^+ \mu^ \mu^+ \mu^) < 1.6 \ (1.2) \times 10^{8}$} and \{${\cal B }(B^0 \rightarrow \mu^+ \mu^ \mu^+ \mu^)< 6.6 \ (5.3) \times 10^{9}$} at 95 % (90 %) confidence level. In addition, limits are set in the context of a supersymmetric model which allows for the $B^0_{(s)}$ meson to decay into a scalar ($S$) and pseudoscalar particle ($P$), where $S$ and $P$ have masses of 2.5 GeV and 214.3 MeV, respectively, both resonances decay into $\mu^+\mu^$. The branching fraction limits for these decays are \{${\cal B }( B^0_{s}\rightarrow SP ) < 1.6 \ (1.2) \times 10^{8}$} and \{${\cal B }( B^0\rightarrow SP )< 6.3 \ (5.1) \times 10^{9}$} at 95% (90%) confidence level.
Feynman diagrams for the $ B ^0_\mathrm{s} \to\mu^+\mu^\mu^+\mu^$ and $ B ^0 \to\mu^+\mu^\mu^+\mu^$ decays, via (a) the resonant $ B ^0_\mathrm{s} \to J/\psi\phi $ SM channel, (b) the nonresonant SM channel and (c) the supersymmetric channel. The latter is propagated by scalar $S$ and pseudoscalar $P$ sgoldstino sfermions. 
SMres.pdf [8 KiB] HiDef png [98 KiB] Thumbnail [124 KiB] *.C file 

SMnonres.pdf [9 KiB] HiDef png [75 KiB] Thumbnail [100 KiB] *.C file 

MSSM.pdf [7 KiB] HiDef png [55 KiB] Thumbnail [111 KiB] *.C file 

(color online). Invariant mass distribution of $K^+ \pi^ \mu^+ \mu^$ candidates. The $ B ^0$ and $\overline{ B }{} ^0_\mathrm{s} $ signal distributions are shown by shortdashed black and longdashed red (gray) lines, respectively. The background shape is denoted by a longdashed black line. The total fit result is shown as a solid blue (black) line. The inset shows the mass distribution centered around the $\overline{ B }{} ^0_\mathrm{s} $ mass. 
jpsiks[..].pdf [94 KiB] HiDef png [289 KiB] Thumbnail [240 KiB] *.C file 

Invariant mass distribution of nonresonant $ B^0_{(s)}\to \mu^+ \mu^ \mu^+ \mu^ $ candidates. The solid (dashed) black lines indicate the boundaries of the $ B ^0_\mathrm{s} $ ( $ B ^0$ ) signal window. The blue curve shows the model used to fit the mass sidebands and extract the expected number of combinatorial background events in the $ B ^0_\mathrm{s} $ and $ B ^0$ signal regions. Only events in the region in which the line is solid have been considered in the fit. 
unblinded.pdf [65 KiB] HiDef png [163 KiB] Thumbnail [147 KiB] *.C file 

Animated gif made out of all figures. 
PAPER2012049.gif Thumbnail 
Combined reconstruction and selection efficiencies of all the decay modes considered in the analysis. The uncertainties shown are statistical. 
Table_1.pdf [46 KiB] HiDef png [169 KiB] Thumbnail [70 KiB] tex code 

Systematic uncertainties on the branching fractions of $ B^0_{(s)}\to \mu^+ \mu^ \mu^+ \mu^ $ . The combined systematic uncertainties are calculated by adding the individual components in quadrature. 
Table_2.pdf [56 KiB] HiDef png [192 KiB] Thumbnail [79 KiB] tex code 
Created on 17 August 2019.Citation count from INSPIRE on 20 August 2019.