The decay of B_s mesons to muon pairs occurs via a "box diagram" in the Standard Model, and is suppressed by the fact that weak interactions must make double work to produce it with respect to more ordinary decay modes. Being weak, weak interactions are unwilling to do that extra work, so the decay probability into muon pairs is of the order of one in a billion.
Several versions of supersymmetric models, on the other hand, would produce a significant enhancement in the production of these fancy decays, whereby the odd-flavoured quark-antiquark pair constituting the B_s meson annihilate by a "neutral current" process, yielding simply a clean pair of muons. The b-quark annihilates with an anti-strange quark more readily in SUSY models because of the possible mediation of extra particles which the Standard Model does not contain.
So, is this a solid indication that Supersymmetry -in some form- is the right theory extending the veteran Standard Model of fundamental interactions ? I doubt it, but I cannot fail to mention that this is indeed an interesting result. I would say that, despite the smaller significance of this result with respect to the one on the muon anomalous magnetic moment (a long-standing discrepancy from the SM which points toward new physics interpretations), it is at the same level of interest. With some caveats.
But let not get ahead of ourselves. What is the result and how does it come about ?
CDF analyzes a total of seven inverse femtobarns of proton-antiproton collisions, sorting out events containing pairs of opposite-sign muons in the range of transverse momentum where a B_s or a B_d signal might contribute. Yes, both B_d and B_s meson decays are studied, because both can yield dimuon pairs, both do so very rarely in the SM, and because the topology is quite similar; moreover, of course both branching fractions are enhanced by new physics intervening in the interaction.
The events are analyzed by a Neural Network which sorts out the most signal-like ones by relying on several observable characteristics of the candidates, notably their probability to originate from a common vertex, the displacement of this vertex from the interaction point, and the kinematics of the muons, plus the isolation of the B candidate they form. Backgrounds to fight against are contributed by combinatorial muon pairs, as well as by B decays to hadron pairs, where the hadrons are misidentified as muons. You can read the details in the paper; here I can only show the main figure of the article, what one might call "money plot". A description follows.
Each row in this 4x3 figure shows a different selection: central-central and central-forward dimuon candidates in the B_s signal search, and central-central and central-forward dimuon candidates summed together in the B_d signal search (central and forward refer to the angle that muons make with respect to the beam line: central means orthogonal to the beam, forward means closer to the beam line; the distinction is only meaningful because the muon detectors in the central and forward regions have different characteristics and performances).
Each row is divided in four histograms of dimuon mass, depending on the value of the neural network output -values closer to 1.0 (the richest in expected signal) being shown in the rightmost panels. The grey histogram shows the expected background level in each selection, the black bullets are the observed counts, and the dark grey histogram is the expected B_s shape from the SM, multiplied by a factor of 5.6 (the best-fit value).
Do not be deceived by the large apparent mismatch between data and background level in the central panels of the first two rows: the event counts are small. CDF chose to not draw error bars (which however with such small event counts mean very little), but you should realize we are talking of few events per bin, so discrepancies are not large. Nevertheless, a small B_s signal might be hinted by the high-NN selections.
CDF computes the p-value of the observed "excess" in the B_s search regions as 0.27% with their standard selection of the NN output; they also quote a 0.67% for a tighter cut which should collect 85% of the B_s signal while getting rid of more background. The latter number is more meaningful to consider in my opinion; CDF however puts 0.27% in the abstract of the paper, so I should oblige, since my name appears in the authors list.
With the observed data, CDF places 95% upper limits on the branching fraction of the B_s and B_d decays to dimuon pairs which are interesting on their own, since they start to affect the region of parameter space where SUSY might contribute. The limits are 6E-9 for the B_d, and only 40E-9 for the B_s -a much looser constraint, due to the existence of an excess of candidates.
If you ask me, the events observed by CDF are a fluctuation, and the branching fraction is indeed what the SM predicts, and not 5.6 times its value (CDF estimates this, but with a large uncertainty: the result for the branching fraction is 18+11-9 per billion). But I am happy to see that something interesting is emerging from these searches. If DZERO and CDF eventually combine their 10/fb datasets on these searches, they might see a signal before LHCb writes a book about it.