Abstract In the framework of the Standard Model we present new theoretical predictions for the branching ratios, double and single differential distributions and forward—backward leptonic asymmetries for the B ¯ s → μ + μ − e + e − decay. In our consideration we take into account the ϕ ( 1020 ) —resonance contribution; the main contributions of four charmonium resonances : ψ ( 3770 ) , ψ ( 4040 ) , ψ ( 4160 ) and ψ ( 4415 ) ; u u ¯ —resonant contribution from ρ ( 770 ) and ω ( 782 ) ; ‘tails’ contributions from J / ψ and ψ ( 2 S ) resonances; non—resonant contribution of the b b ¯ - pairs, bremsstrahlung and the contribution of the weak annihilation. We provide the prediction for the branching ratio of B ¯ s → μ + μ − e + e − decay with and without the ϕ ( 1020 ) —resonance contribution. We use the model of vector meson dominance (VMD) for calculation of resonances contributions and take into account all substantive terms in B ¯ s → μ + μ − e + e − amplitude that was not considered in the previously papers.