  The rare decays $B_{s,d}^0 \to \ell^+\ell^-$ are analyzed within the general
framework of the aligned two-Higgs doublet model. We present a complete
one-loop calculation of the relevant short-distance Wilson coefficients, giving
a detailed technical summary of our results and comparing them with previous
calculations performed in particular limits or approximations. We investigate
the impact of various model parameters on the branching ratios and study the
phenomenological constraints imposed by present data.
