Motivated by the first LHCb searches for the rare \( {\overline{B}}_{s,d}^0 \) → J/ψμ+μ− decays, we perform a detailed study of these processes within the QCD factorization formalism. Since the transverse size of the J/ψ meson is small in the heavy quark mass limit, this formalism is generally expected to hold for these decays. We include both the leading- and next-to-leading-order QCD corrections to the hard-scattering kernels, which are convoluted with the light-cone distribution amplitudes (LCDAs) of the initial- and final-state hadrons. It is numerically found that, depending on the model parameters for the leading-twist B-meson LCDA, the maximum branching ratios of \( {\overline{B}}_s^0 \) → J/ψμ+μ− and \( {\overline{B}}_d^0 \) → J/ψμ+μ−, integrated over the dimuon invariant mass squared q2 from 1 GeV2 to ( \( {m}_{B_{s,d}} \) − mJ/ψ)2, can reach up to 2.21 × 10−9 and 7.69 × 10−11 at the leading order in αs, respectively. After incorporating the non-factorizable one-loop vertex corrections, these branching ratios are further reduced by about one order of magnitude, with \( \mathcal{B}{\left.\left({\overline{B}}_s^0\to J/\psi {\mu}^{+}{\mu}^{-}\right)\right|}_{q^2\ge 1{\textrm{GeV}}^2} \) = 2.88 × 10−10 and \( \mathcal{B}{\left.\left({\overline{B}}_d^0\to J/\psi {\mu}^{+}{\mu}^{-}\right)\right|}_{q^2\ge 1{\textrm{GeV}}^2} \) = 1.07 × 10−11. In addition, we have presented the dimuon invariant mass distributions of the individual and total helicity amplitudes squared, as well as the differential and integrated longitudinal polarization fractions of the J/ψ meson, which could be probed by the future LHCb and Belle II experiments with more accumulated data.
